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Executive Summary 


NASA/ Lewis Research Center is sponsoring a program for providing computer codes for 
analyzing and designing turbomachinery seals for future aerospace and engine systems. The 
program is made up of three principal components: 1) the development of advanced 3-D 
Computational Fluid dynamics codes 2) The production of simpler 2-D industrial codes and 3) 
the development of a Knowledge Based System(KBS) that contains an expat system to assist in 
seal selection and design. 

The 3-D code is being produced by a major subcontractor. Computational Fluid Dynamics 
Research Corporation (CFDRC) of Huntsville ,AL., who are enhancing an existing CFD code, 
REFLEQS. The first task of CFDRC has been to concentrate on cylindrical geometries with 
straight, tapered and stepped bores. Improvements have been made by adoption of a colocated 
grid formulation, incorporation of higher-order, time-accurate schemes for transient analysis 
and high-order discretization schemes for spatial derivatives. This report describes the 
mathematical formulations and presents a variety of 2-D results, including labyrinth and brush 
seal flows. Extensions to 3-D are presently in progress. 

Three industrial codes have been produced which are capable of being run on a PC. 

• SPIRALG predicts performance characteristics of gas-lubricated, spiral-groove 
cylindrical and face seals including eccentricity and misalignment ( four degrees of 
freedom which consist of two orthogonal displacements and two orthogonal angular 
misalignments), which represent extensions to the present state of the art. The code 
produces seal loads and moments, minimum film thickness, axial flow, power loss 
and up to thirty two frequency dependent cross-coupled spring and damping 
coefficients. Arbitrary placement of grooving and the dam region is permitted as well 
as user selection of the spiral-groove pumping direction. The code is coupled to an 
optimization code that will allow for determination of optimum groove geometry on 
the bases of stiffness , pumping capacity and flow. A code option is the use of a 
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Romberg extrapolation procedure for producing rapid and accurate results. Included in 
this report is a comprehensive theoretical development of the code, several sample 
cases and validation information. 

• The industrial code ICYL is intended for use in analyzing cylindrical seals operating 
with incompressible fluids. The code includes film turbulence, and inertia effects at 
inlet and exit and at boundaries where sharp clearance discontinuities result, such as 
hydrostatic recesses. Configurations include plain circular, hydrostatic, multi-lobe, 
tapered and Rayleigh Step geometries. An important feature of the code is the 
incorporation of roughness on the seal housing or rotating shaft. It thus permits 
analysis of damping seals which are finding favor in advanced cryogenic 
turbomachines. The code produces seal loads, and righting moments, flows , power 
loss, clearance and pressure distributions , up to thirty two cross-coupled dynamic 
spring and damping coefficients as well as critical mass and frequency. This report 
describes the theoretical development , includes examples of code usage and 
validation against other codes and information in the literature. 

• The industrial code GCYL analyzes cylindrical gas seal configurations. 
Configurations include plain circular, hydrostatic, multi-lobe, tapered and Rayleigh 
Step geometries. The code produces seal loads, and righting moments, flows , power 
loss, clearance and pressure distributions , and up to thirty two frequency dependent, 
cross-coupled dynamic spring and damping coefficients. This report describes the 
theoretical development , includes examples of code usage and validation against 
other codes and information in the literature. 

The functions of the Knowledge Based system are 1) to integrate the scientific and industrial 
codes 2) to provide a user friendly graphical user interface and 3) to include an expert system for 
seal selection, analysis and design. A significant requirement is portability between a PC and 
UNIX based workstation. The two operating systems selected are OS/2 with the Presentation 
Manager interface for the PC and UNIX with OSF/MOTIF for a workstation. A two track 
development effort is being pursued. The scientific codes are being developed under UNIX and 



the industrial codes are being developed using the OS/2 operating system. The user interfaces are 
being developed using object oriented tools and C++ which are portable between OS/2 and 
UNIX. The initial development platform will be OS/2 and porting to UNIX will be 
accomplished by recompilation. This report discusses development plans and presents some 
OS/2 graphical user interfaces accomplished with the industrial codes. 

In addition to code and interface development, the project requires technology transfer to both 
government and non-government facilities. A peer panel has been established, whose function is 
to guide program development, and annual workshops are held to transfer information. The first 
workshop was held on March 26, 1991. 
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1.0 INTRODUCTION 


NASA's advanced engine programs are aimed at progressively higher efficiencies, greater 
reliability, and longer life. Turbomachinery for future aerospace engine systems will require 
advanced seal configurations to control leakage, control lubricant and coolant flow, prevent 
entrance of contamination, inhibit the mixture of incompatible fluids, and assist in the control of 
rotor response. 

A seven year program has been devised with the objective of providing to NASA and the U.S. 
Aerospace Industry, three dimensional scientific codes and simpler industrial codes for analyzing 
and designing optimized advanced seals with minimal development time. 

The program provides three interdependent parallel paths: 

1. The development of scientific Computational Fluid Dynamics (CFD) codes capable of 
producing full three-dimensional flow field information to enhance understanding of 
flow phenomena and mechanisms, to contribute design guidance for complex situations, 
and to furnish accuracy standards for less sophisticated analyses. All tasks involving 
three-dimensional code development will be accomplished by a major subcontractor to 
MTI, CFD Research Corporation (CFDRC) 

2. The development of industrial codes for expeditious analysis, design and optimization of 
turbomachinery seals. The industrial codes will consist of a series of separate, 
stand-alone codes that will be integrated by a Knowledge Based System (KBS). 

3. The development of expert systems to assist users to select an appropriate seal type for 
their application, provide design guidance, and assist in interpreting data from the 
analysis programs. 

The analysis codes and the expert systems developed by the three activities will be integrated 
into a unified system by the KBS which will provide access to and link all the various 
components. The key features of the KBS include the following: 
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• Access to all analysis codes and expert systems 

• An easy- to- use, consistent user interface for all KBS components 

• Utility functions such as printing and browsing output files 

• Plotting of output data from analysis programs 

• Database of analytical models and other supporting information 

• Portability between PC and Workstation environments 

An important aspect of the contract is technology transfer to the industrial, government and 
academic communities. This is being accomplished through annual workshops, reports, and code 
distribution through NASA. The first workshop was convened on March 26, 1991 with over 65 
attendees. A Peer Review Panel has also been established consisting of seal experts and 
cognizant representatives from industry, government and the academia. The Peer Panel provides 
technical guidance to the program. 

This report covers the effort completed during the first year of the program which included the 
following: 

• The development of advanced algorithms and validation of the CFD codes with 
emphasis on cylindrical geometries. 

• Delivery of three industrial codes to NASA for Beta testing 

• A gas lubricated spiral groove code SPIRALG for analyzing spiral groove 
cylindrical and face seals. 

• A cylindrical incompressible seal code ICYL for analyzing a wide variety of 
cylindrical geometries including roughened surface seals. 

• A cylindrical compressible seal code GCYL for analyzing a wide variety of 
cylindrical gas seal geometries. 

• The establishment of a detailed plan for the implementation of the Knowledge Based 
System using the PC and OS/2 as the principal up front interface and operating system 
respectively. Several industrial codes were implemented and information is presented 
in this report. As a result of the first workshop, a need for a UNIX operating system 
was expressed by the attendees and the members of the Peer Panel. Resolution is 
presently being accomplished by KBS software that will be portable to both operating 
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systems and will be operable utilizing PCs or Workstations. The PC can only act as 
an interface for the scientific codes that must reside on a mainframe or workstation. 
The industrial codes can be self contained in a PC environment. The universal 
approach to the KBS will be presented at the next workshop and descnbed in the next 

annual report. 
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2. 0 COMPUTATIONAL FLUID-DYNAMICS DEVELOPMENTS 


2.1. Introduction 

2.1.1 Ppyelopment of the cy lindrical seals CFP code 

The objective of Task I is to develop a three-dimensional CFD code for analysis of 
flows in straight, tapered and stepped cylindrical seal configurations. This code will 
be capable of solving three-dimensional Navier-Stokes equations in generalized, 
body-fitted coordinates with provisions for polar and cylindrical systems. The 
features which are relevant to the seals program include. 

1. Stationary and rotating coordinate systems; 

2. Steady-state and time-accurate solution capability; 

3. Advanced turbulence models for high shear rotating flows; 

4. Incompressible and compressible flow solutions; 

5. Variable physical properties (viscosity, density, specific heat , efc.); 

6. Cavitation effects; 

7. Provision for stepped surfaces and injection ports; 

8. Inclusion of viscous dissipation and phase changes in energy equation; 

9. Treatment of sources due to external fields, e.g. electromagnetic and 
electrostatic; 

10. Variable surface roughness treatment; 

11. Provision for effects of pre-swirl and upstream effects; and 

12. Customized input and output features for cylindrical seals. 

The code will utilize solution procedures and schemes that are accurate, efficient 
and robust to include all these characteristics for high-aspect ratio computational 
cells typically encountered in seed geometries. 
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2.1.2 Focus of W ork for the First Yeai 

During the past year the focus of work has been to develop, implement and test 
several new concepts in the basic code REFLEQS. The two-dimensional version of 
REFLEQS was selected as the starting point for all the development work. Two- 
dimensional problems are sufficiently general, so that once proven, the concepts can 
be extended to three dim ensions in a straightforward manner; at the same time the 
complexity of the coding is sufficiently low so that rapid development and incorpo- 
ration of these concepts in the 2-D code are possible. Following are the modifica- 
tions and improvements which were made in the basic REFLEQS code. 

1. Adaptation of a colocated grid formulation in which the velocity compo- 
nents as well as the scalars are stored at the computational cell center as 
against the earlier staggered grid formulation where the velocity compo- 
nents are stored on the cell faces; 

2. Use of Cartesian velocity components as the primary velocity variables in 
place of the velocity projections which were used before; 

3. Incorporation of high-order time-accurate schemes for transient flow 
analyses which include a) PISO algorithm, b) Crank-Nicholson method, 
and c) three-point second-order backward time-differencing method; and 

4. High-order discretization schemes for spatial derivatives. These include, 
in addition to central differencing, third-order upwind-biased scheme, 
Osher-Chakravarty scheme, and minimod limiter scheme. 

Some of the items described above merit further attention at this point, and the 
merits of these and the reasons for implementation are discussed below. 

2.1.2.1 Colocated grid formulation w ith Cartesian Components. Figure 2.1 illustrates 
the velocity-pressure location arrangement for staggered and colocated grid configu- 
rations. 
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Figure 2.1 Examples of Staggered and Colocated Grids 

A number of CFD codes currently in use are based on the staggered grid approach. 
These codes use finite volume methods with segregated or coupled equation 
solution methods. The chief reason for the use of staggered grid approach is to 
avoid the phenomenon of the odd-even decoupling of pressure when solving 
■incompressible flows. By locating the velocity nodes on the cell faces instead of the 
cell centers, these velocities now can be linked directly to the pressures at the two 
nearest cell-centers. This provides a strong coupling between the velocities and the 
pressures and avoids the checkerboard pressure pattern. In recent years, however, 
interest in colocated grid formulation has been renewed 1 ' 3 . Coupled with this 
approach is also the use of Cartesian components as the primary velocity compo- 
nents. This combination offers a number of advantages which are listed below. 

1. A common control volume for mass, momentum and energy conserva- 
tion eliminates many calculations which are repeated for the various 
control volumes in the staggered grid, e.g. evaluation of the link coeffi- 
cients that are needed to set up the discretized form of the flow equations 
except the mass equation; 

2. In a single grid cell, the number of interpolations required to calculate the 
velocity components at the cell faces is minimized; 

3. In a staggered grid, the boundary condition implementation is more 
involved due to the physically displaced control volumes. In a single grid 
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cell the boundaries are the same for all variables, and the treatment is 
uniform; 

4. The use of Cartesian components ensures that the treatment of body-fitted 
coordinate (BFC) and complex grids is relatively simple and easy to 
understand and implement. Use of Cartesian components also simplifies 
the procedure of velocity interpolation to cell faces. Use of any other 
primary velocity variable can make the interpolation procedure very 
cumbersome in a complex geometry, since the local angle relations must 
be taken into account; and 

5. Implementation of higher-order spatial discretization schemes is simpler 
with the single cell approach. 

The main drawback of the colocated grid formulation is that the coupling between 
the pressure and velocities cannot be maintained as easily as in the staggered 
approach. Recently, however, several methods have been proposed and used 
successfully to overcome the problem of odd-even decoupling of pressure. The 
particular formulation used in the present work is discussed in detail in a later 

section. 

717 ? piso Algorithm. The earlier formulation in REFLEQS for transient flow 
calculations involved several iterations of the overall solution procedure for each 
time step. The solution procedure thus could become expensive since, in effect, 
each time-step solution involved the solution of the corresponding steady-state 
solution. The PISO algorithm is designed to calculate transient flows with a non- 
iterative scheme. The algorithm consists of a predictor step where an intermediate 
solution is calculated, followed by a series (typically 2 or 3) of corrector steps which 
improve the accuracy of the predicted solution. Several of the steps in the overall 
scheme are implicit, so that the algorithm is much more stable with respect to the 
time step size as compared to an explicit time-marching scheme. Due to the non- 
iterative nature of the PISO algorithm the overall computational costs for this 
method can be substantially smaller than the iterative methods. The algorithm can 


2-4 



also be set up to achieve higher order time-accuracy. For these reasons, the devel- 
opment of PISO algorithm for transient analyses was considered. 

2.2 Mathematical Formulation 

In this section a discussion of the theoretical details of the various procedures 
implemented in the 2-D code is given. The basic differential and finite difference 
equations for the fluid flow are shown. A discussion of the mass interpolation 
procedure used in the code follows next, and finally the solution steps needed in the 
two basic algorithms: SIMPLEC and PISO are given. 


2.2.1 Flow Equat ions and Discretization 

In the Cartesian tensor form the fluid flow equations can be written as 


Continuity 


momentum 


d(P«i) _ 0 
at dxi 

. dXj 




a_ 

dXi 


dx, 


+ Bf 


or, alternatively, the general transport equation for any flow property, <t> f is 


( 2 . 1 ) 

( 2 . 2 ) 



(2.3) 


where B,- is the body force, T is the diffusion coefficient, and S* is the source term 
associated with the variable <p. This source term, then, would contain the pressure 
terms and other body force terms for the momentum equations. 

The flow equations are next transformed to a generalized, Body-Fitted Coordinate 
(BFC) system which allows the grid to conform to the problem geometry. The 
switch to the BFC system (^.q) is done using the transformation 
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«-«(x,y) 

ri=ri(x,y) 


(2.4) 


The partial derivatives can then be transformed as: 


8 , a 8 

_=^— + 77 x — 

8y 8rj 


8,8 8 


(2.5) 


The two-dimensional general transport equation for <p then becomes 


a ( p/») t a(pu») , a(py») _ » 

St 3| 3n 3$ 


;r(v|-v{)^J + ^/r(vi7-V7))|i] + /(s' + s*) 


( 2 . 6 ) 


where U and V are the contravariant velocity components and S' is the source term 
associated with curvilinear nonorthogonal part of the viscous stress tensor. 


U = /(&cU + S y v) 


(2.7) 


V=/( TJx« + 


s' = ±[,r( v 4.vn)g] + A 


/r(v7iv$)± 


( 2 . 8 ) 


(2.9) 


with u and v as the Cartesian velocity components along x and y directions. S ' is the 
additional source term which is generated during the transformation of the 
diffusion terms; it is zero for an orthogonal grid. / is the Jacobian of the transforma- 
tion. 


d(s,y) _ 

afeu) " 


x s 


yri 


( 2 . 10 ) 
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The transformed continuity equation is 


a(/p) . a(pLr) a(pv) _ 0 ail) 

dt k an 

For the solution methods used in the present work the transport equations are 
integrated over a general computational cell in the grid to generate algebraic 
equations which link the variables in the cell with those in the surrounding cells. 



Figure 2.2 Computational Control Volume Grid and Nomenclature 
The discretized transport equation for a variable <p then is given by 

(«* - p $ pv is ♦ s* 


( 2 . 12 ) 


where y is the volume of the cell .the subscript nb denotes all neighboring cells and 
a n bf etc. are the link coefficients which consist of the convective and diffusive terms 
linking the cell centers with those around it. and the form of the coefficients 
depends on the method used for spatial differencing. Thus, e.g., for the upwind 

scheme. 
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a E<p = + max M pU)e Sri e ) 


(2.13) 


apf is obtained by summing all the link coefficients of neighboring cells. 

The transport equation. Equation 2.12 is written for u and v velocities and solved 
sequentially to update the values of the components to u* and v*. 

The updated velocity components at this point do not satisfy the continuity 
equation. To impose this the flow variables are assumed to have a form 


u n+l -U* + U 

(2.14) 

v n+l _ + v ‘ 

(2.15) 

pit+1 + p' 

(2.16) 

(f 1 * 1 sp‘ + p' 

(2.17) 


The momentum equations are rewritten for the velocity corrections as 


| fl p„ + “j u p = X a nbu u nb " (**«£ V$ + ^ui} Vi] 

[ a Pv + V P = S, a nbv v'nb * v\ Zy + *vi\ Vi] Vy) 


(2.18) 

(2.19) 


where the subscript for pressure denotes partial derivatives. The method of 
treatment of the term under the summation sign decides the algorithm which is being 

used. If u nb = 0 is assumed the summation term is simply dropped and a SIMPLEC 
type algorithm results. In the SIMPLEC algorithm the individual corrections at 

neighbors are taken to be the same as at point P, i.e. u nb - u'p- With this approxima- 
tion the summation term is merged with the term on the left hand side. In the 
PISO algorithm, discussed in section 2.2.3, the neighbor velocity corrections u nb are 
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evaluated with the last iteration (correction) level of the predictor-correctors 
procedure. 


The continuity equation is then integrated over the typical computational cell to 
give 


£_ + (pile) A? sin a e - (pU w ) A w sin a w + ( pV J), A„ sin a n - (pV \ A s sin a s m 

At 


where 


vt — (pli)e Ag sin otg ~ (plf)ai A ip sin Ota; + (pV)n An sin otn * (pl^)s A$ sin 0% 


v n 



Figure 2.3 Control Volume Nomenclature for BFC Grids 


( 2 . 20 ) 


Referring to Figure 2.3, a denotes the angle between the constant £ and Tj lines, and 
the subscripts w,e,n,s refer to the sides of the cell, and A denotes the area of a cell 
side. The correction in the contravariant components are 


Li’ = u'Zx + v €y ^ 2 ’ 21 ^ 

v' = u‘t} x + v rjy (2.22) 
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Expressions for U ' and v'. Equations 2.21 and 2.22 are now substituted in Equation 
2.20 to provide a pressure-correction equation of the form 


f 

“7 + a Pp ?P ~ 5- Pnb ' m 

At 

(2.23) 

For incompressible flows the p' term is taken as zero. 

For compressible flows it is 

expressed as 


• V P • 

(2.24) 

P RT p P 

and is absorbed in the coefficient associated with p’ p . 

Solution of Equation 2.23 

provides pressure corrections at the cell centers which then are used to calculate 

corrections in other variables: 

« = " P% Zx + dut] Prj %) 

(2.25) 

v =- [dv$ P$ + Pn ty) 

(2.26) 

• P 

(2.27) 

P ~ RT 

Finally these corrections are vised to update the velocities, pressure and the density. 

u " +1 = u* + u 

(2.28) 

V n+1 sb v* + V 

(2.29) 

pn+1 _ p« + p 

(2.30) 

pn+l - + p‘ 

(2.31) 


This completes a typical iteration in the SIMPLEC procedure for steady-state 
equations. For steady solutions the time-step At provides one form of underre- 
laxation. The PISO algorithm also follows similar steps, and at this point the 
predictor and the first corrector step in this algorithm would be complete. 
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Evaluation of the mass source term in Equation 2.20 is critical to the success of the 
colocated grid formulation. To calculate the mass source term the contravariant 
components at cell faces are needed which in turn are interpolated from the cell- 
centers. Improper interpolation procedures can lead to odd-even decoupling of 
pressure. The formulation used in the present work is described in the following 
subsection. 

2J2J2 Mass-Carryi ng Velocity Interpolation 

The basic concepts of this interpolation procedure will be developed for a 1-D 
problem for ease of understanding. The procedure outlined is similar to that in 
Reference 2. Extension to two dimensions is described next. 



41 1S/1 sms 


Figure 2.4 Schematic Grid for Mass-Carrying Velocity Interpolation 

For this 1-D problem, the contravariant component, U, is the same as the Cartesian 
component u. The value of u at the cell face 'e' is to be calculated using the veloci- 
ties at P and E. Use of a simple average to get: 

«e = ^(« P + M E ) < 2 -3 2 > 

leads to odd-even decoupling of pressure. This is due to the fact that with this 
definition, the velocities at the cell faces do not directly depend on the pressure 
difference between the neighboring cell centers. This strong coupling is achieved in 
the staggered grid by physically shifting the location of cell face e. In the colocated 
grid formulation, the cell face velocities must be calculated with a direct coupling 
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with neighboring pressures. To achieve this, the momentum equations for the 
velocities at P, E and e are written: 


“ p =(X<h,“Jp-'*p(i;) J( +s p 

(2.33) 

“ £ =(X^“J e -^(^) e +s e 

(2.34) 

“, = (X a 'rtu S « 

(2.35) 


where a nb = ^ / 4 = ^ , etc. 

[Here the prime denotes that the term is divided by the coefficient associated with 
the velocity at P or E or e.] At this point the effect of time step and/ or underre- 
laxation is ignored; it will be added at a later stage. 

The various terms in Equation 2.35 are then approximated using the corresponding 
terms in Equations 2.33 and 2.34. Thus 

(X a nbu u nb + s )e “ 2 [£ + S ) P + E + S ) £ ] (236) 

and 

d; = l(d p + d £ ) (2.37) 

Substitution in Equation 2.35 yields 

u e = \ [(X a 'nbu u nb + S)p + (X a nbu u nb + S )e] ’ 2 E ) 

Expansion of the terms under the summation signs gives rise to a large connectivity 
with nodal velocities and is complicated to evaluate. Instead, these terms are 
expressed as 
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(2.39) 


a nbu u nb + 
(S a nbu u nb + 


s)p - Up- dp j- 

s) £ = u E -d l 


M 

Bxjp 

M 

dxk 


(2.40) 


Using these expressions in Equation 2.38 one obtains 


u « ~ 2 [ Mp + 



which is much simpler to evaluate. If a further approximation is made: 


(2.41) 


1 



/dp\ + /3p| 

2 

p |0x)p E )e. 

2 

\dxjp [dx/E, 


(2.42) 


the final form of u e is obtained: 



(2.43) 


The other cell face velocity, u w now can be calculated in a similar fashion. When 
these velocities are used to calculate the mass source term, the pressure derivatives 
add together to generate a fourth-order derivative of pressure. This serves to 
suppress odd-even decoupling of pressure by providing a stronger coupling between 
cell pressures. 


Next step is to include the effect of time derivatives and/or underrelaxation terms. 
Improperly done, this can give rise to steady-state solutions which depend on the 
size of the underrelaxation used. The following analysis is similar to that developed 
in Reference 4. Equations 2.39 to 2.41 are rewritten as 


a Pu + 


pV 

At 




(2.44) 
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(2.45) 



(2.46) 


where it is assumed that the time term serves as the underrelaxation factor for 
steady solution. These equations now have a modified coefficient with the veloci- 
ties, and also include a term with last iteration /time-step values of the velocities, 
denoted by the superscript o. 


Following an analysis similar to that given earlier, the expression for u e is now 
written as 


- 2 [ + “E ] ' ic IS ' 2 IS + IS] 


(M 


At 


a Pu + 


pV' 
At > 




(2.47) 


and the term associated with the time term is calculated as 



(2.48) 


Equation 2.48 is the final form of the interpolation procedure that can be used to 
calculate cell face velocities. Extension of these concepts to two-dimensional BFC 
grids is discussed below. 

A typical computational cell in a 2-D BFC grid is shown in Figure 2.5. The contra- 
variant components, U and V along the % and tj axes are to be calculated using the 
velocities at the surrounding cell centers. 
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Figure 2.5 Convective Fluxes Passing Through BFC Control Volume Cell Faces 


The components are given by 

U e = (uJt x + vtfy)e (2.49) 


V s = {uJVx + u/Tj y )s 


(2.50) 


Thus, to arrive at the cell face contravariant components, both u and v have to be 
interpolated at each face. Consider the component at face e . To calculate the 
Cartesian components the momentum equations at cell centers P and E are used. 
These are 



(2.51) 


(2.52) 


where terms such as now contain the metric coefficients which multiply the 
pressure derivatives. Similar expressions for these components at cell center P can 
be written. Since \l e is the component along the $ direction, the pressure redistribu- 
tion is applied only in the £ direction. Thus, 
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(2.53) 


l(— ) - 1 

[dp] + / 3 p\ 

2 . 

.We WpJj 


(pv). 


At 


l pu 


pv 1 

At > 


«?4( u 2 + M $) 



(Cp + v E ) - d v ^ e 


11— \ - 1 

73p\ + idp\ 11 

ml 2 

.U/E \d£)p\ 


(pvt 

dt 



l_(ro.l( r | + *$) 


df ^ 


(2-54) 


Finally, the contravariant component U e is calculated using Equation 2.49. The 
component along the rj direction at, e.g., the south face is evaluated using the 
Cartesians at cell centers P and S, and the pressure redistribution terms along T) 
direction only are taken into account. 


This is the basic formulation which has been incorporated in the 2-D code to avoid 
the odd-even decoupling of pressure in the colocated grid formulation. Extension of 
this concept to a 3-D grid is straightforward, and not discussed here. 


2.23 Solution A lgorithms 

A description of the solution algorithms currently implemented in the 2-D code is 
given in this subsection. The major solution steps are outlined for four schemes: 1) 
Modified SIMPLEC, 2) Crank-Nicholson, 3) Three-point backward time-differencing, 
and 4) PISO. The first three schemes use essentially the same basic solution algo- 
rithm, while the PISO scheme requires additio n al steps. 

1. Modified SIMPLEC: This algorithm consists of three main steps. 

a. Evaluation of an intermediate velocity field, u* and v* by solving 
momentum equations such as 2.12 with lagged pressure terms; 
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b. The continuity equation. Equation 2.20 is solved next for pressure 

corrections. These are used to update pressures and velocities. Density is a 
updated for compressible flows; 

In the standard SIMPLEC procedure this marks the end of an iteration and step (a) is 
taken next. In the modified SIMPLEC procedure, a step is added to ensure a tighter 
continuity condition. 

c With the updated flow variables in step (b) the mass source term in 
Equation 2.20 is reevaluated. Using the new mass source term the continu- 
ity equation is solved for pressure corrections. These are then used to 
update the flow variables again. The link coefficients in the pressure 
equation are kept frozen during this step. Step (c) is repeated till a suitable 
criterion is reached. At this point one iteration is considered complete 
and the next iteration started with step (a). 

Steps (a) to (c) are repeated till a suitable convergence criterion is reached. This 
algorithm is the default option for steady-state flow solutions. 

2. Crank-Nicholson Scheme: This algorithm is adopted for transient flow analysis 
and is formally second-order accurate in time. This is achieved by evaluating all 
convective and diffusive fluxes at time level (n + 1/2) where n is the old time level. 

The algorithm consists of the following steps: 

a. Evaluate all flux terms using the last time-step variable values, i.e. at level 
n. These fluxes are not updated till the next time-step is taken; 

b. Intermediate velocity field u* f v* is calculated using Equation 2.12 with 
lagged pressure. The convective and diffusive fluxes are calculated using 
the following expression. 

Ifc + M = « Ifc + f D t * I 1 - «>{fc + ! d T < 2S5 > 
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where the superscript k denotes the iteration level, a is called the Crank- 
Nicholson parameter, and controls the implicitness of the scheme as well 
as the time-accuracy. The scheme is second order in time for a = 0.5;. 

c Pressure corrections are evaluated using Equation 2.20. The mass source 
term in this equation is evaluated using 

m ■ am * + (1 - a)m n (2.56) 

If needed, iterations on this step are done by updating the mass-source 
term. This is done in a fashion outlined in step (c) of the SIMPLEC 
procedure; 

d. Steps (b) to (c) are repeated till a convergence criterion is reached. At this 
point the solutions at iteration level k are taken to be solutions at new 
time level n+1; 

e. Time is advanced by a step, and calculations are started at step (a). 

By changing the value of a, the order of accuracy and nature of the scheme can be 
changed. Thus a = 1 corresponds to the Euler backward time discretization which is 
implicit in time and first-order accurate. 

3. Three-point backward time-discretization: This is another second-order time- 
accurate method for transient flow analysis. The high-order accuracy is achieved by 
discretization of the time derivative using a three-point method. 

3ft _ 3ft W+1 - 4ft w + ft”' 1 (2.57) 

& lAt 

where superscripts n-1, n, and n+1 denote different time levels. This is also an 
iterative algorithm and the steps are: 
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a. Time terms are calculated based on the expression given above; 

b. This step essentially involves iterations of the SIMPLEC procedure with 
the additional time term at n-1 carried as a source; and 

c At convergence of step (b), the variables are updated and a new time step 
is taken. Calculations for this step start from step (a). 

4. PTSO Algorithm; This is the non-iterative algorithm implemented in the 2-D 
code for transient flow analysis. The solution steps consist of a predictor step 
followed by a series of corrector steps; the basic procedure is oudined in a paper by 
Issa 5 . The algorithms for incompressible and compressible flows differ somewhat 
and each is outlined below. 


Incompressible flows. For these flows there is no density variation, so that the 
energy equation need not be calculated during the predictor-corrector sequence for 
velocities and pressure. The algorithm steps are: 


a. With a new time-step, an intermediate velocity field, u* and v* is 
calculated using momentum equations. Equation 2.12. This is the 
predictor step; 

b. For the first corrector, the pressure correction equation. Equation 2.20 is 
solved to yield pressure corrections. These are then used to update flow 
variables to u**, and p*; 

c A pressure correction equation followed by an explicit momentum 
corrector equation completes the second corrector step. The pressure 
correction equation is 


a Pp Pp~^ a nbp Pnb m " X a nbu ( u nb ~ w nb) “ 


V" 1 / ++ * \ • * 

Z, a nbv\°nb' v nbl' P=P ~P 


(2.58) 


The pressure corrections are used in the momentum correction equation 
for w' 
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*Pu 


. v 1 / ** * \ . ' dp] • »»• ** 

W P ~ \ u nb ~ u nbl ~ ^P hj^j ' u p~ u p ~ u p 


(2.59) 


where the tilde denotes that the time term has been absorbed in that term. 
A similar equation for v' is solved. With these corrections, the flow 
variables are updated to u***, v***, and p**. This completes one time step 
in this algorithm. The corrected variable values are taken as the new time- 
level values; and 

d. A new time step in taken and calculations started at step (a). 

Compressible Flows. PISO algorithm for compressible flows is more involved, since 
the density variations have to be calculated. This is incorporated in the momentum 
and continuity equations. In addition, the energy equation has to be solved at each 
corrector stage to update the temperature. The solution steps in this method are: 

a. With a new time-step the first predictor and corrector steps are taken 
using the procedures outlined in Equations 2.12 and 2.20. The values of 
the variables at this stage are u**, v**, p*, and p*; 

b. Using the updated values, the energy equation is assembled and solved to 
generate the corrected values of the temperature T*; 

c The second corrector step now involves solution of p' again using a 
correction equation. The form of this equation is: 


a Pp Pp - X a nbp Pnb " * m "Xp a nbu(u„i,- U nb )- 

At 

X P* a nbv (v*l - V* b ) + b, p' = p”-p* 


( 2 . 60 ) 


Solution of this equation, the p' field is then used to update the velocity 
fields to «*** and v*** using equations such as: 


~ *** v ** i(dp I ^ 

a Pu U P ~ 2* a nbu u nb " “ l"0^~ I + 


(2.61) 
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A similar equation for v*** is used. At the end of this step, the flow 
variables are updated to u***, v***, p** and p**; 

d. The energy equation is again assembled and solved to yield the corrected 
values of the temperature, T*\ This is the second energy corrector step; 

e. This completes the so-called two stage scheme. If desired the steps (d) and 
(e) can be repeated to add more corrector stages. At the end of this series, 
the updated values of the flow variables are taken to be the new time level 
values; and 

f. A new time step is taken and the calculations are started at step (a). 

2.3 Status of the 2-P Code 

In the present form, the 2-D code is based on the colocated grid formulation de- 
scribed above, and uses the Cartesian components as the primary velocity variables. 
The modified SIMPLEC algorithm is the default for steady-state flow solutions, 
while the schemes available for transient analyses are: 

1. First-order accurate backward differencing; 

2. Three-point, second-order accurate time differencing; 

3. Crank-Nicholson scheme; and 

4. PISO scheme. 

A comprehensive set of boundary conditions is provided which includes: 

1. Specified velocities; 

2. Specified pressure; 

3. Wall boundary; 

4. Symmetry condition; and 

5. Zero-gradient extrapolation condition. 
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Profiles for all the flow variables can also be specified at all boundaries if known, e.g. 
for pre-swirl inlet conditions in a seal flow. 


In the colocated formulation, the pressure at cell faces is used to calculate the 
pressure derivatives needed in the momentum equations. Proper evaluation of 
pressures at the boundaries thus becomes important. Two types of second-order 
accurate pressure extrapolation procedures are used at the boundaries in the 2-D 
code. 



Figure 2.6 Nomenclature for Pressure Boundary Condition Interpolation 

1. For an inflow condition the slope at the boundary is assumed to be the same 
as at the nearest cell center. Referring to Figure 2.6, the boundary pressure can 
be calculated as 

(2 Ax p + 2Ax c) Ax p 

Pwh-*- FPp’l E — rP E 

(Ax p + Ax p ) (ax p + Ax £ ) E (2.62) 

2. For extrapolation, symmetry, and wall conditions, the slope at the boundary is 
calculated using a second order extrapolation. The slope at the boundary is 
given by: 


dp _p p Ax%- p<; Ax* - p n ( 2Ax p + Ax - Ax$ 


(2Ax p + Ax s f Ax p - Ax j (lAx p + Ax s ) 

(2.63) 
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If the slope is specified to be zero, the expression for the pressure at the boundary 
becomes 


Pn = 


Pp + Axj 

(2Ax p + Ax s f - Ax} 


(2.64) 


For swirling flows the centrifugal acceleration has to be taken into account. This is 
accomplished by specifying the pressure slope in terms of the angular speed co at the 
boundary. 


|?) =(py a \ (2.65) 

\oyln 


Pn 


p p (2Ar p + Ax s f - p s t (pyo^ln dx p (2Ar p + Ax s ) 


(lAx p + Ax s f -Ax} 


3 AXp + Ax s 


( 2 . 66 ) 


Finally, if a body-force is present, such as gravity, it is used to specify the pressure 
derivative. Thus, e.g. for gravity force, the pressure derivative and boundary 
pressure are given by 



= -(Pg)n 


(2.67) 


Pn 


_ p p (2Ar p + Ax s f - Ax} (pg)n Ax p (2Ax p + 


( 2 . 68 ) 


(lAxp + - Ax} 


3 Ax p + Ax s 


The code is capable of handling incompressible and compressible flow. Several 
turbulence models are incorporated which are 1) mixing length model, 2) low 
Reynolds number k-e model, 3) standard k-e model with wall functions, and 4) 
multi-scale k-e model. 


One of the problems specific to seal geometries is a rotor undergoing whirling 
motion in a seal as shown in Figure 2.7. The 2-D code can be used to simulate such a 
rotor with the assumption that the axial pressure gradient is zero, or in other words, 
when there is no leakage. To facilitate this, a coordinate frame whirling with the 
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rotor is selected. When the whirl orbit is circular, this transformation reduces the 
time-dependent problem to a quasi-steady one. The momentum equations are then 
solved in terms of the relative velocities. Rotation of the axes gives rise to addition- 
al source terms, the Coriolis and centrifugal accelerations, which are added to the 
momentum equations. 


+ (convective terms ) « (Pressure + diffusive terms) - pQ x(fl x ;)-2pHx« 

dt 

Centrifugal Coriolis 


(2.69) 


where Q is the whirling angular speed. Finally, the velocity boundary conditions at 
the rotor and stator wall have to be modified, and are given by 


stator: u = - cox r 

rotor: (f2- ©) x r-[{2- cojx A 

where £ is the position vector joining the centers of the rotor and the stator, and a> 
is the angular speed of the rotor. This formulation also has been incorporated in the 
2-D code, and can be invoked by specifying a non-zero angular speed for precession. 


I 



Figure 2.7 Schematics of Rotor /Stator Configuration With Circular Whirling Orbit 
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2.4 2-D Code T est Results 


The 2-D code has already been used to calculate a range of standard flow solutions. 
These computational test cases were designed to assess overall accuracy of the code 
as well as the accuracy of the various physical models. Presented in this section are 
solutions for a number of selected test cases which have a direct relevance to the 
seals application. Accuracy of the numerical results is checked against analytical 
solutions in several cases. Two seal calculations are also presented at the end of the 
section which are checked against experimental data. These test cases serve to prove 
the capability of the computer code to calculate accurate and physically sound 
solutions. 
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2.4.1 Flow in an Annulus Between Two Cylinders 
Problem Specification 

• Developing flow in an annulus between two cylinders. 

• Narrow annulus, ratio of inner to outer radii = 0.995. 

• Laminar flow, Reynolds number based on outer radius = 100. 

• Slug flow at inlet, fully-developed at the exit. 

Benchmark Data 

• Analytical solution for fully-developed flow. 

Grid 


• 20 cells in both axial and radial directions, evenly spaced. 

• A maximum aspect ratio of 3xl0 4 . 

Boundary Conditions 

• Uniform flow at the inlet. 

• Fixed pressure at the outlet. 

• Wall conditions at both cylinder surfaces. 

Results 

• Flow description is given in Figure 2.8a. 

• Figure 2.8b shows the calculated axial velocity as a function of radius at the 
last axial station. The analytical solution for the fully-developed flow is 
also plotted for comparison. Excellent agreement between numerical and 
analytical results is obtained. 
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(a) Flow Geometry 



(b) u Velocity 


Figure 2.8 Flow in an Annulus Between Two Cylinders. Rj /Rq = 0.995, 20 x 20 grid, 

grid aspect ratio = 3x 104, Re^ = 100 
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2AJZ Flow Between Rotating Cylinders 
Problem specification 

• Row in the annulus between two cylinders. 

• Inner cylinder rotating at 28650 rpm., stationary outer cylinder. 

• Narrow annulus, ratio of inner to outer radii = 0.995. 

• Laminar flow, no flow in axial and radial directions. 

Benchmark data 

• Analytical solution for the stable Taylor-Couette flow. 

Grid 

• 4 cells in the axial and 50 cells in the radial direction. 

• Maximum aspect ratio = 3.6xl0 4 . 

Boundary conditions 

• Periodicity in axial direction, solutions at first and last axial stations are 
taken as identical. 

• Wall conditions with specified angular speed =28650 rpm. at inner 
cylinder. 

• Wall conditions with zero angular speed at outer cylinder. 

Results 

• Row conditions and geometry shown in Figure 2.9a. 

• Computed tangential velocity as a function of radius shown in Figure 2.9b. 
Corresponding analytical solution also plotted for comparison. 
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(b) Tangential Velocity w 


Figure 2.9 Flow Between Rotating Cylinders, Ri/Ro = 0.995, cOq = 0, C0i = 28650 rpm. 

4 x 50 grid, aspect ratio = 3.6 x 10 4 
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2.4.3 Couette Flow 
Problem specification 

• Planar flow between two parallel, infinite plates. 

• Top plate moving at a constant speed. 

• Uniform pressure gradient is applied. 

Benchmark data 

• Analytical solution for Couette flow. 

Grid 


• 3 cells in the flow direction, 20 cells across the gap in the plates with even 
spacing in both directions. 

Boundary conditions 

• Periodicity conditions imposed at the cross-planes. 

• Stationary wall at bottom plate. 

• Wall condition with specified velocity at top plate. 

Physical models 

• The pressure gradient term is included as a special source term in the 
main-flow momentum equation. 

Results 


• Flow geometry and parameters are shown in Figure 2.10a. 

• Flow solutions are obtained at several values of pressure gradient parame- 
ter ranging from -3 to +3. The corresponding numerical and analytical u 
velocity profiles are shown in Figure 2.10b. 


2-30 



1 I 

-0 3 0 0 0 5 


(b) u Velocity as a Function of P 


Figure 2.10 Couette Flow. Pressure Gradient Coefficient P = 


2|iU \ 9x 
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2.4.4 Planar Wedge Flow 
Problem specification 

• Laminar flow in a narrow wedge-shaped passage. The top block is at rest; 
the bottom plate is moving. 

• Flow passage is very narrow (L/h = 3X10 3 ) 

Benchmark data 

• Analytical solution for the planar wedge flow. 

Grid 


• BFC grid with 192 cells along the length and 40 cells across the gap, evenly 
spaced. 

Boundary conditions 

• Wall condition on the slider block. 

• Wall condition with specified velocity on the bottom plate. 

• Specified pressures at the two passage openings. 

Results 

• Grid in shown in Figure 2.11a. 

• Streamline pattern in the passage in shown in Figure 2.11b. 

• Comparison of computed and analytical u velocity profiles at several 
locations along the length are shown in Figure 2.11c 

• Pressure across the passage is constant; computed and analytical pressure 
profiles along the length are shown in Figure 2. lid. 
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(a) 192 x 40 Computational Grid 
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(b) Streamlines 


Figure 2.11 Planar Wedge Flow. Length = 0.1m, Height: Minimum = 3 x 10' 5 m, 
Maximum = 3 x lCHm, y Scale Enlarged 200 times in (a) and (b) 
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2.4.5 Laminar Flow Over Backstep 

Problem specification 

• Laminar flow over a planar backward facing step; expansion ratio = 1:1.94. 

• Reynolds number = 100, 300, 389, 500 and 648. 

Benchmark data 

• Experimental data of Armaly, et al. 6 . 

Grid 

• 110 cells in the flow- wise direction, 40 cells across. Cells clustered near the 
step in flow direction, and in the passage upstream of the step. 

Boundary conditions 

• Specified uniform axial velocity at inlet; value of the axial velocity varied 
depending of the Reynolds number. 

• Specified pressure at outflow boundary. 

• No-slip at all wall boundaries. 

Results 

• Reattachment length for the flow as a function of the Reynolds number is 
plotted in Figure 2.12a. Computed results are compared with experimen- 
tal results 6 . 

• Figure 2.12b shows the streamline pattern for Reynolds number = 500. 
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(a) Reattachment Length L = x/h 
as a Function of Reh 



(b) Streamlines, Reh = 500 
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Figure 2.12 Laminar Backstep Flow. 110x40 Grid 



2.4.6 Turbulent Flow in a Plane Channel 
Problem specification 


• Turbulent flow in a planar channel; Reynolds number = 61,600. 

• Treated as a developing flow problem, with fully developed flow at the 
channel end. 

Benchmark data 

• Hot-wire measurements by Laufer 7 . 

Grid 


• 50 cells in the flow direction with even spacing. 40 cells in the cross 
direction with clustering near the wall for a specified cell width. 

Boundary conditions 

• Uniform flow specified at inlet. 

• Constant pressure at the outflow. 

• Wall conditions at upper and lower walls. 

Physical models 

• Standard k-e model for turbulence with wall functions. 

Results 

• Flow details are shown in Figure 2.13a. 

• Computed profiles of turbulence kinetic energy and streamwise velocity at 
the last station are plotted in Figures 2.13b and 2.13c. Also shown in these 
figures are the experimental data from Laufer 7 . 
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(b) Turbulent Kinetic Energy ( c ) Mean Velocity 


Figure 2.13 Turbulent Flow in a Plane Channel. 
Re = UoH/ y = 61,6000. Data from Laufer^ 
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2.4.7 Turbulent Flow Induced by Rotating Disk in a Cavity 
Problem specification 

• Calculation of the flow induced by a rotating disk in an enclosed cavity. 
Benchmark data 

• Experimental measurements from Daily and Nece 8 . 

Grid 

• 40 cells in the axial direction, 60 cells in the radial direction with clustering 
near the walls. 

Boundary conditions 

• Specified angular velocity for the rotor walls. 

• Wall conditions for all other boundaries. 

Numerics and physical models 

• Central differencing with 0.05 damping. 

• Standard k-e model with wall functions. 

Results 

• Flow geometry as shown in Figure 2.14a. 

• Normalized radial and tangential velocities at a given radius are shown in 
Figures 2.14b and 2.14c. Also shown in the figures are the experimental 
data from Daily and Nece 8 . 
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(a) Flow Geometry 



(b) Radial Velocity (c) Tangential Velocity 


Figure 2.14 Turbulent Flow Due to a Rotor in an Enclosed Cavity. 
Experimental Data from Daily and Nece 8 
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2.4.8 Flow Between Stator a nd Whirling Rotor 
Problem specification 


• Flow in the clearance between a stator and a whirling rotor. 

• Circular whirl orbit assumed. Calculations are performed in a coordinate 
frame whirling with the rotor. 

• Solutions computed at whirl speeds of 0.01, 0.5 and 1 times the shaft 
angular speed. 

Benchmark data 

• None. 


Grid 


• 40 cells in the circumferential direction and 10 cells in the clearance, 
evenly spaced. 

Boundary conditions 

• Wall conditions with wall velocities corresponding to the transformed 
frame. 

• Cyclic conditions assumed in the circumferential direction. 

Numerics and physical models 

• Central differencing with 0.05 damping. 

• Standard k-e model for turbulence. 

Results 

• Grid and flow geometry shown in Figure 2.15a. 

• Pressure distribution in the clearance shown at the three whirl frequencies 
in Figures 2.15b through 2.15d. 
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(a) 40 x 10 Computational Grid 
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(b) Pressure Contours for Synchronous Whirl, Cl = co 


Figure 2.15 Flow in the Seal Clearance for a Whirling Rotor. 
Stator Radius = 50.1478 x lO^m, Clearance = 495pm, e = 0.8 x Clearance, 
Shaft Rotation Speed to = 5085r|m, Clearance Exaggerated for Clarity 
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(c) Pressure Contours, Subsynchronous Whirl ft = 0.5© 
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(d) Pressure Contours, Subsynchronous Whirl ft = 0.01© 


Figure 2.15 Flow in the Seal Clearance for a Whirling Rotor Continued 
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2.4.9 flow Qver a Bank of Tubes (Bru sh Seals) 
Problem specification 


• Planar flow over a bank of tubes. This flow is similar to that in a brush 
seal. 

Benchmark data 

• None. 


Grid 

• Three rows of tubes with three tubes in each row considered. 

• 60 cells in both directions; a BFC grid is employed. 

Boundary conditions 

• Specified uniform velocity at the inlet. 

• Specified pressure at the outflow. 

• Symmetry conditions specified at the two remaining outer boundaries. 

• Wall conditions specified on all tube surfaces. 

• Tubes simulated using blocked cells. 

Results 

• Grid and flow geometry is shown in Figure 2.16a. 

• Computed pressure contours for this flow are shown in Figure 2.16b. 
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2.4.10 Annular Seal Flow 
Problem specification 

* Calculation of turbulent flow in an annular seal. 
Experimental data 

• Experimental data by Morrison, et al. 9 . 

Grid 


• 25 cells in the radial direction, 58 cells in the axial direction; cells in radial 
direction clustered near the walls. 

Boundary conditions 

• Experimental profiles of the velocities and turbulence quantities at inlet 
boundary. 

• Specified pressure at the outflow boundary. 

• Wall condition with specified angular speed at rotor wall. 

• Stationary wall conditions at stator wall. 

Numerics and physical models 

• Central differencing with 0.01 damping. 

• Standard two equation k-e model for turbulence. 

Results 

• Geometry of the rotor is shown in Figure 2.17a, and the experimental 
setup is shown in Figure 2.17b. 

• Computed and experimental contours of the axial, azimuthal and radial 
velocities are shown in Figures 2.18, 2.19 and 2.20, respectively. 

• Figure 2.21 shows the computed turbulent kinetic energy profiles. 
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•4 JO mm 
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(a) Geometry of Annular Rotor 



(b) Annular Seal Inlet Geometry 


Figure 2.17 Flow Details for the Annular Seal. Reynolds Number Based on the 
gap = 27000, Taylor Number = 6600, Shaft Speed = 3600 rpm 
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Experimental 


Numerical 


Figure 2.18 Annular Seal Flow. Contours of the Scaled Axial Velocity u x /U. (a) 
Near Inlet (0 < x/c < 4), (b) Near Outlet (25 < x/c < 29.4) 


{at —u —.3.1 —0.1 —xi —0.1 —c.1 —0.1 — ai —ai — 


Experimental 


Numerical 


Figure 2.19 Annular Seal Flow. Contours of the Scaled Azimuthal Velocity , 
Ue/wshaft (a) Near Inlet (0 < x/c < 4), (b) Near Outlet (25 < x/c < 29.4) 










Figure 2.20 Annular Seal Flow. Contours of Scaled Radial Velocity u r /U. 

Near Inlet Details (0 < x/c < 4) 




(a) 


<b) 


Figure 2.21 Annular Seal Flow. Contours of Scaled Turbulent Kinetic Energy, 
l-(uj 2 + Uq + u x 2 )/U 2 - Numerical Results, (a) Near Inlet (0 < x/c < 4), 

^ (b) Near Outlet (25 < x/c < 29.4) 
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2-4.11 Seven Cavity Labyrinth Seal 
Problem specification 


• Calculation of turbulent flow in a seven-cavity labyrinth seal. 

Experimental data 

• Experimental data by Morrison, et a/. 10 . 

Grid 

• 30 cells in the axial and radial directions per cavity. 

• 10 cells in the radial clearance between the rotor tooth and the stator. 

• Stretching used to cluster the grid near the rotor and stator walls. 

Boundary conditions 

• Experimental profiles for velocities and turbulence quantities at inlet 
boundary. 

• Specified pressure at outflow boundary. 

• Wall condition with specified angular velocity at rotor walls. 

• Wall conditions at stator wall. 


Numerics and physical models 

• Central differencing with 0.01 damping. 

• Standard two equation k-e model for turbulence. 

Results 

• Details of the rotor are shown in Figure 2.22a, and the experimental setup 
is shown in Figure 2.22b. 

• Computed and numerical velocity vector plots are shown in Figure 2.23. 

• Computed and experimental contours of the axial, radial and tangential 
velocities are shown in Figures 2.24, 2.25 and 2.26, respectively. 

• Figure 2.27 shows computed contours of the turbulent kinetic energy. 
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(a) Details of the Labyrinth Rotor 



(b) Seal Inlet Geometry 


Figure 2.22 Seven Cavity Labyrinth Seal Flow. Reynolds Number Based on the 
Clearance = 28000, Taylor Number = 7000, Shaft Speed = 3600 rpm 
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Experimental 


Numerical 


Third Cavity 

Figure 2.23 Seven Cavity Labyrinth Seal. Velocity Vector Plot 
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Experimental 


Numerical 




Third Cavity 


Figure 2.24 Seven Cavity Labyrinth Seal. Contours of Scaled Axial Velocity, u x /U 
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Experimental Numerical 

Third Cavity 

Figure 2.25 Seven Cavity Labyrinth Seal. Contours of Scaled Radial Velocity, u r /U 
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Experimental Numerical 


Third Cavity 


Figure 2.26 Seven Cavity Labyrinth Seal. 
Contours of Scaled Azimuthal Velocities, u 0 /w S h a ft 
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First Cavity Third Cavity 


Figure 2.27 Seven Cavity Labyrinth Seal Flow. Contours of Normalized Turbulent 
Kinetic Energy l( u ;2 + U( 2 + u ;2) /u2 . Numerical Resulls 


2-56 



REFERENCES 


1. Peric, M., Kesseler, R., and Sheurer, G., "Comparison of Finite Volume Numeri- 
cal Methods with Staggered and Colocated Grids/' Computers and Fluids, Vol. 
16, No. 4, 198 8 , pp. 389-403. 

Z Rhie, C. M., and Chou, W. L., "Numerical Study of the Turbulent Flow Past an 
Airfoil with Trailing Edge Separation," A1AA /., Vol. 21, 1983, p.1525. 

3. Parameswaran, S., and Sun, R., "Numerical Simulation of Turbulent Flows 
Around a Car-Like Body Using the Non Staggered Grid System," AIAA Paper 89- 
1884, AIAA 20th Fluid Dynamics, Plasma Dynamics and Lasers Conference, 
Buffalo, NY, June 12-14, 1989. 

4. Majumdar, S., "Role of Underrelaxation in Momentum Interpolation for 
Calculation of Flow with Nonstaggered Grids," Numerical Heat Transfer, Vol. 
13, 1988, pp. 125-132. 

5. Issa, R. I., "Solution of the Implicitly Discretised Fluid Flow Equations by 
Operator-Splitting," J. of Comp. Phys., Vol. 64, 1985, pp. 40-65. 

6. Armaly, B. F., Durst, F., Pareira, J. C. F., and Schonung, B., "Experimental and 
Theoretical Investigation of Backward-Facing Step Flow," J. Fluid Mech., Vol. 
127, 1983, pp.473-496. 

7. Laufer, J., "Investigation of Turbulent Flow in a Two-Dimensional Channel," 
NACA Report 1053, 1951. 

8. Daily, J. W., and Nece, R. E., ' Chamber Dimension Effects on Induced Flow and 
Frictional Resistance of Enclosed Rotating Disks," Trans, of ASME, J. of Basic 
Eng., Vol. 82, pp.217-232. 

9. Morrison, G. L., Johnson, M. C., and Tatterson, G. B., " 3-D Laser Anemometer 
Measurements in an Annular Seal," ASME Paper 88-GT-64, ASME Gas Turbine 
and Aeroengine Congress, Amsterdam, The Netherlands, June 6-9, 1988. 

10. Morrison, G. L., Johnson, M. C., and Tatterson, G. B., " 3-D Laser Anemometer 
Measurements in a Labyrinth Seal," ASME Paper 88-GT-63, ASME Gas Turbine 
and Aeroengine Congress, Amsterdam, The Netherlands, June 6-9, 1988. 


2-57 


3.0 Industrial Code SPIRALG - Gas-Lubricated, Spiral Groove, Cylindrical and Face 
Seals 

Spiral groove bearings and seals are used to provide stability, load support and pumping for both cylindrical 
and face seal geometries. In the case of a cylindrical seal, grooves are usually designed to pump against 
each other In a symmetric arrangement to provide enhanced stabBity. A lightly loaded cylindrical seal 
operating at a low compressibility number will produce a force that is nearly 90 degrees out of phase with 
the displacement which will tend to destabilize the rotating shaft The Introduction of spiral grooves can 
significantly increase the component of force In phase with the displacement and decrease the out of phase 
component thereby improving stability. 

In the case of a face seal or thrust bearing, spiral grooves are often Introduced as the primary means of load 
support Since a symmetric arrangement Is not possible In a radial geometry, the grooves are usually 
designed to pump towards an ungrooved dam region. The resistance of the dam region increases as the 
film thickness decreases hence the pumping pressure rise Increases thereby giving rise to a positive axial 
stiffness. The spiral grooves can also be used to pump against an applied pressure gradient thereby 
resulting in either reduced or reversed leakage. 

The computer code SPIRALG predicts performance characteristics of gas lubricated, spiral-groove, 
cylindrical and face seals. Performance characteristics include load capacity, leakage flow, power 
requirements and dynamic characteristics in the form of stiffness and damping coefficients in 4 degrees of 
freedom for cylindrical seals and 3 degrees of freedom for face seals. These performance characteristics 
are computed as functions of seal and groove geometry, loads or fBm thicknesses, running speed, fluid 
viscosity, and boundary pressures. 

The basic assumptions that have gone Into the computer code are listed below: 

1. The flow is assumed to be laminar and isothermal. 

2. Inertial effects are neglected. 

3. The gas is assumed to be ideal. 

4. The film thickness is assumed to be small compared with seal lengths and diameters but large 
compared with surface roughness and the mean free path of the gas. 
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5. Narrow groove theory is used which characterizes the effects of grooves by a global pressure 
distribution without requiring computations on a groove by groove basis. This involves neglecting 
edge effects and local compressibility effects associated with groove to groove pressure variations. 
In genera), narrow groove theory is valid when there are a sufficiently large number of grooves so 
that 2rslns/N g < < 1, where p is the groove angle and N g Is the number of grooves. 

6. Transient effects are treated with the use of small perturbations on a primary steady state flow. 
These transient effects are characterized by stiffness and damping coefficients that are dependent 
on the disturbance frequencies. 

7. Although displacements and misalignments are treated, machined surfaces for face seals are 
assumed to be flat and machined clearances for cylindrical seals are assumed to be constant. 

The above assumptions still leave the code applicable to a broad range of applications. Seals generally 
have small clearances and gasses have low densities resulting In sufficiently low Reynolds numbers for 
laminar flow. Practical designs should contain a fairly large number of grooves to ensure smooth, Isotropic 
operation. At high sealed pressure differences, the flow could become sonic thereby invalidating the first 
two assumptions but this will usually not be the case and can readily be checked based on the predicted 
leakage flow. Elastic and thermal distortions as well as machining tolerances should also be estimated to 
validate the constant clearance assumption. The overall accuracy of the program will depend on the grid 
size used. Factors such as high compressibility or squeeze numbers, small values of the minimum film 
thickness to clearance ratio and large values of the length to diameter ratio could require either a large 
number of grid points or carefully selected variable grids. 
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3.1 Theoretical Development 


The first formulation of the equations governing gas lubricated spiral groove bearings is generally credited 
to Vohr and Pan [1 ]*. A more concise formulation is given in a second report [2] by these authors that has 
been used by Smalley [3] as a starting point in his generalized numerical treatment of the performance of 
spiral groove gas bearings. The work performed by Smalley may be applied to both bearings and seals. 
A principal limitation in all of the above references relates to the fact that solutions have only been provided 
for one dimensional forms of the equations which have been obtained by linearizing them based on near 
concentric and aligned conditions. The work described here deals with the numerical solution of the 
nonlinear equations for gas lubricated spiral groove seals at both eccentric and misaligned conditions. 


Formulation of equations governing gas lubricated spiral groove seals 


For completeness, a derivation of the narrow groove equations for spiral groove gas bearings and seals 
along the lines of that developed in Reference 2 will be provided here. Coordinate variables will be used 
to make the equations applicable to both cylindrical and face seals as can be seen with the aid of Figure 
3-1. The circumferential coordinate, 6, is as shown In Figure 3-1. The transverse coordinate Is described 
by the variable, s, which is taken to equal the radial coordinate, r, for a face seal and the axial coordinate, 
z, for a cylindrical seal. The quantity r, when It appears wDI denote radial position for a face seal and should 
be set equal to the shaft radius, Rg for a cylindrical seal. 


The isothermal, compressible form of the "Reynolds* equation may be written as a flow balance equating 
the divergence of the flow vector, qf , to the flow per unit area squeezed out by the time rate of decrease 
of the film thickness, 


tf-q' 



r 80 


Qa 


(3-1) 


The local flow vector q' = q^1 + q'J represents the mass flow rate per unit transverse length divided 
by the density at a reference pressure, which may be written in vector form as 


* Numbers in brackets refer to references given at the end of this section. 
""Nomenclature is given at the end of this section. 
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Figure 3-1 Coordinate system for spiral groove analysis. 
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Since surface motion will be in the circumferential direction, the surface velocity vectors may be written as 
u, = nw.,1 and u 2 = r» 2 l and the components of q' become 
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(3-2) 
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(3-3) 


The "squeeze film" term or displaced mass flow per unit area due to film motion, divided by the density at 
Po is 

qi - - (3-4) 

One could substitute Equations (3-2) - (3-4) for the corresponding flow quantities in Equation (3-1) to obtain 
the usual form of the compressible Reynolds Equation which could In principle be solved, for any film 
thickness profile, h(s,6) and appropriate boundary conditions, for the pressures or flow components to 
obtain the pressure distribution. These could in turn be integrated to obtain the various forces and moments 
associated with the given bearing geometry. The torque opposing the motion of say the smooth surface 
may be determined, once the pressure distribution is known, by integrating the shear stress relationship that 
arises in the development of Reynolds equation 
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The difficulty encountered in obtaining full numerical solutions to the above equations relates to the 
complexity of the grid network necessary to adequately describe the geometry of a surface containing the 
large number of spiral grooves usually required to provide sufficiently smooth pressure distributions to make 
the load characteristics independent of whether shaft displacement Is over a ridge or over a groove. Narrow 
groove theory is generally used to circumvent this difficulty (References 1 - 3). It will be implemented here, 
as well and Is described below. 


Narrow groove theory provides the limiting form of the solution to Equations (3-1) - (3-5) as the number of 
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Figure 3-2 Schematic of spiral groove parameters, global and local pressures. 
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grooves, N Q , becomes large, with the groove angle, p and the groove to pitch ratio, a , held constant. The 
discontinuities In film thickness associated with the grooves wUI give rise to discontinuities in the pressure 
gradients at the ridge-groove interfaces as illustrated schematically in Figure 3-2. The local pressure profile 
p is shown by the sawtooth lines whose lower vertices, for purposes of illustration, are connected by the 
"global* pressure profile, p. The global pressure profile does not necessarily lie at the lower vertices of the 
local pressure profile but could lie anywhere between the lower and upper vertices. In the limit as the 
number of grooves becomes large the curve connecting the upper vertices will approach the curve 
connecting the lower vertices. This limiting behavior is not true of dp /BO, dp fit or h, which will have 
different values over lands and grooves no matter how large the number of grooves. Narrow groove theory 
requires the development of expressions for the local (primed) quantities in terms of global quantities that 
approach single valued limits as the number of grooves becomes large. The local fDm thickness and 
pressure derivatives over the grooves will be denoted by h g , dip /d0) g and (3p'/3s) g respectively and 
by h r (9 p /0Q) r and (9p'/9s) r over the ridges. (The subscript r has been used here to denote ridges 
for consistency with References 1 - 3 and should not be confused when used in a different context later to 
denote the right hand boundary pressure or with the radial position variable, r, which is not used as a 
subscript.) 

When the number of grooves becomes large, the sawtooth portion of the localpressure variation may be 
approximated with linear representations as shown in Figure 3-2. Thus, equating pressures over a groove- 
ridge pair in the circumferential direction 


AE . «. *£l . fel + M 

A0 A6 A0 ldeJ g A6 { 00 ), AO ’ 

Noting that AO^AO = a and A0 f /A0 ■ 1 - a and replacing Ap/A0 with dp/dO as A0 - 0, the above 
equation becomes 


ie . +( i -*)(&) . 

30 laoj. 1 ; la0j r 


(3-6) 


The corresponding relationship in the transverse direction, 


f -*(£)/<’-*>(£), • 


(3-7) 


is obtained in a similar manner. 
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The remaining two equations required to solve for the four local pressure derivatives are obtained from 
continuity considerations. 

First, the pressure must be continuous at each groove-ridge interface, thus the derivative of the pressure 
In the direction of the interface, Vp'»tp, must also be continuous. The second requirement is for 
continuity of the flow normal to each groove-ridge interface as measured in a frame of reference moving with 
the grooves, (q' - ru 1 hp r /p 0 i)»ry. The unit tangent and normal vectors for a logarithmic spiral are 
given by 


\ = costfl + sirys] , rip = sirysl - costfj . 


The first of the above conditions requires continuity of 


oosii^ ♦ shpie: 
r ae as 


at each groove-ridge interface or 

The second condition requires continuity of 


(q^-ru>h-£^)slnp - q'cosp . 
Po 


One may substitute Equations (3-2) and (3-3) for the circumferential and transverse components of the flow 
vector at each groove-ridge interface, respectively to obtain 

(3-9) 

The density variation term, p'/p 0 is continuous at each interface and cancels out of Equation (3-9). 
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Equations (3-6) - (3-9) represent the four linear equations needed to solve for the local pressure derivatives. 
We may obtain the solution by first solving Equations (3-8) and (3-9) for the components of the local 
pressure gradient over the grooves in terms of those over the ridges. The resulting equations may be 
written as 



hgCO**p + h 3 sin 2 p 1 1 h 3 - h, 3 
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(3-11) 

- 6|tr(<i> 2 -unstop cos p . 

h # 

One may now substitute Equation (3-10) for (3p'/30) f in Equation (3-6) and Equation (3-11) for 
(dp'/ ck)^ in Equation (3-7) to obtain 2 linear equations for the components of the ridge pressure gradient 
which may in turn be solved to yield the following expressions: 
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The components of the local groove pressure gradient may ge expressed in terms of the above ridge 
components by simple rearrangement of Equations (3-6) and (3-7): 


l(ap ; ] . 1 - a 1 [ap'] + 
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(3-14) 
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fe) = -1.- c (j^) + liE . (3-15) 

v as J a « [ as ), a as 

Now that expressions have been developed for the components of the local pressure gradients in terms of 
global ones, it is necessary to determine the global flow components q # and and a global squeeze film 
term q A that may be substituted for the local ones in the flow balance given by Equation (3-1). These global 
flow components are determined by matching mass flow rates over a groove-ridge pair with the mass flows 
obtained by Integration of the local flow components over the same interval. 

If 0 g Is taken as the circumferential coordinate at the start of a groove, the transverse flow crossing an arc 
at fixed s, subtending a groove-ridge pair in the interval 0 g < 0 < 0 g +A6 is given by the left hand term 
in the relationship 


«,-Ae 

/ q,'rde 




12 n p 0 



h* p ( dp'j 
12|i Pplae j r 


rA6 r 


q,rA0 . 


The approximation to the integral in the above expression was obtained by dividing the Integration Interval, 
AO into sub-intervals for the groove, A0 g and ridge, A0 r and approximating q' a , noting that as the number 
of grooves becomes large 3p'/3s, will approach a constant value within each sub-interval. Since the 
pressure at the groove-ridge interface is continuous, the local density variation term, p' /p 0 was replaced 
. by its global value p/p ff The far right hand term in the above expression is based on the definition of the 
transverse component of the global flow rate described above. The right two equalities may be solved for 
q, as 


q. 


12 p p 0 lasj. 


hf* p 
12»i Pa 


(1 " a) fe) * 

V“/r 


( 3 - 16 ) 


One may obtain a relationship for the circumferential flow component in a similar manner by Integrating 
O'®, given by Equation (3-3), at fixed 0, over a groove-ridge pair ( s g < s < s g +As ), approximating the 
integral over each sub-interval as above and equating the result to c^A s. The resulting expression may be 
written as 


Qe “ 


hg 




12|» p 0 r(aeL 12p p 0 


h? 


f(1-«) 


r v 50 Jr 


+ r ^«JE. loh (1 . a)hrJ # (3.17) 


Po 


By integrating the squeeze film term q A ', given by Equation (3-4), over an area rABAs, equating it to 
q A rA0As and noting that the groove area fraction will be a and the ridge area fraction will be 1 - a, the 
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following expression is obtained: 


q A “ -°>M) * (3-18) 

The global shear stress t, may be determined by integrating the local shear stress t', given by Equation 
(3-5), with respect to 6 over the Interval 0 g < 0 < 0 g +A0, invoking the narrow groove approximations and 
equating the result to tA0. The resulting expression may be written In the form 


T 



♦ 


2 r 



+ |ir(« 2 



(3-19) 


Equation (3-1) may now be applied directly to the global flow vector q = c^l + q # j, as V«q = q A and 
by substituting Equation (3-18) for q A and putting the result in dimensionless form one obtains: 


tf-Q 


1 d 
RdS 


(RQ.) 


, 130 . 
R dO 


-i-l(«**H r )(1*P)] . 
ot 


(3-20) 


The components of the global flow vector, and q # are given In terms of the local pressure derivatives 
by Equations (3-16) and (3-17) respectively. These local derivatives are. In turn, given in terms of the global 
ones by Equations (3-12) - (3-15). The global flow components may be expressed completely in terms of 
global pressure derivatives by first substituting Equations (3-14) and (3-15) for the local pressure derivatives 
over the grooves and then substituting Equations (3-12) and (3-13) for the local pressure derivatives over 
the ridges. One may then collect terms and put the resulting two equations In dimensionless form to obtain 
the following expressions for the components of the dimensionless flow vector Q = c^l + Qj: 


Q. «-(1*P) 


44 * H) 


P ♦ Ajk^Rsinp - A(«8 + H r )R 


(3-21) 


Q. 


-<1*P) 



♦ la h) p " A * k < ncosp • 


The dimensionless variables associated with the above equations are 


(3-22) 


The dimensionless gage pressure P in the above equations is taken relative to the absolute pressure p 0 
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p.£a. a.jast*. n-A. s.i-.R.f.i.i.r.i .<« 3 > 

Po C s Po c Rg Rq 2A h, 

which will henceforth be taken as the minimum of the two boundary pressures in absolute units. The 
dimensionless parameters associated with the above equations are 


GiiuRq 


I* 


■ , A 4 ■ A6 wa(1 - e)sinp , S> - — — — , 5 » , a * — (3-24) 

PoC 2 « C l r ♦ l 8 


and the column matrix containing spiral groove coefficients, kfofi F), in the above equations is 


a(1 -aMl^-DWp+I* 

(1 - «)r® ♦ c 

a(1 - aXT 8 - losing cosp 
(i - «)r* + « 

«(i - g)(r* - i^oos 2 ^ +r* 

(i - «)r* ♦ « 

Slz l 1 -■ 

(1 - «)r* + a 

r 

(T-I)slnp 

(1 

a(1 - «)(!* - 1)(T - 1)slnp cosp 
(1 - a)r* ♦ a 

a(1 - alfl* - 1)(T - POOS 2 ? ♦ af ♦ (1 - a)! 32 
(1 - a)r* ♦ a 


(3-25) 


Only the first 4 components of k are used in Equations (3-21) - (3-22). The remaining components are used 
in evaluating the shear stress. The relationships for k,, is 1,2, 3, 4 derived here are consistent with 
Equation (3.27) of Reference 2. 


The global shear stress is obtained by substituting Equation (3-14) for pp'/36) 0 /r in Equation (3-19) then 
substituting Equation (3-12) for (9p'/30) r /r in the resulting expression. The latter result may be expressed 
in dimensionless form as 
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The equations presented thus far are directly applicable to either a cylindrical seal or a face seal. As 
mentioned earlier, a face seal is represented in the above equations by setting the transverse coordinate s 
equal to the radial coordinate r. This is equivalent to setting S * R In dimensionless form. A cylindrical seal 
is represented in dimensionless form by setting S = Z and R>1. 

The quantities required to characterize the groove dimensions are shown In Figure 2. If by convention u 
is taken to be positive (surface motion in the direction of increasing 6), then the groove angle, p , will be 
the angle measured from the groove to the direction of surface motion associated with o. A positive acute 
value of/> will tend to pump in the positive S direction If the grooves are on the stator and in the negative 
S direction for grooves on the rotor. By setting the groove depth parameter 4 = 0, r becomes 1 and 
Equations (3-20) - (3-26) reduce to those for ungrooved seals. By treating k and 4 as sectionally continuous 
functions of S, these equations may be applied to composite smooth and grooved geometries with Q # and 
P held continuous at all transition boundaries. 


The film thickness relationship for H r which may be applied to either a cylindrical seal or a face seal Is 


H r - 1 - e, - (€ x + tS)cos0 - (ey -4>S)sin6 (3-27) 

with e z a 0 for a cylindrical seal and c » c « 0 for a face seal. 

The boundary pressures will be taken to be p, and p r at the Inside and outside radii respectively for a face 
seal or at the two ends (z = -L/2 and Z = L/2) for a cylindrical seal. This is expressed in dimensionless 
form 

P - P, at S - S, , P - P r at S - S r . (3-28) 

The remaining boundary condition relates to periodicity with respect to 8 which requires P and to have 
the same values at 8 * o as they do at 8 = 2*: 


p le-o “ PI.- 2 , «nd Q,|,^ - 0,1,.*. . (3-29) 

The above treatment is intended to represent a complete statement of the mathematical problem for 
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determining the pressures and surface shear stresses in plain or spiral groove face or cylindrical seals. The 
rest of this section will deal with the numerical determination of the pressure distribution and the 
computation of related quantities such as loads, leakage, power loss, stiffness and damping. 

Discretization of pressure equations 

Discretization will be carried out with the use of the ceil method (4] which involves the performance of a flow 
balance about each interior grid point One may integrate Equation (3-20) over an arbitrary control area 
‘within a seal 


/ tf-QdA + / 4l(««*H r )(1 + P)]dA = 0 

A A 

and apply the divergence theorem to the first integral on the left to obtain the relationship 


i Q-ndS + f +H.)(1 +P)JdA - 0 , (3 " 30) 

which will be used as a starting point In the discretization process. 

A grid network may be set up along with flow control areas about each grid point as shown in Figure 3-3. 
The grid will contain M lines In the S direction including boundaries and N lines in the 0 direction from 
0 = 0 to 0 = 2k , inclusive. The grid points at the Intersections of these lines are noted by the solid circles. 
Row control areas to be used in evaluation of the Integrals In Equation (3-30) are set up about each grid 
point as shown by the shaded area In Figure 3-3. The comers of the flow control area denoted by the 
shaded points marked 1,2,3, and 4 are located at the geometric centers of the rectangles formed by the grid 
lines and will be referred to as half grid points. The flow components labeled Q j 2 etc., represent the 
components of the flow vector in the positive coordinate directions as indicated by the arrows. The 
subscripts (12 etc.) refer to the line connecting points 1 and 2. and the superscripts (+,-) refer to the 
positive or negative side of the point of Intersection with the grid line. 

We will adopt the convention that the subscripts i,J refer to grid points and subscripts such as I+V4.1+V4 refer 
to half grid points. The value of the radius R at half grid point 2 would thus be R^. The differential control 
length, dS, in Equation (3-30) will be approximated by the lengths of the various lines or arcs bounding the 
flow control area thus ASj 2 refers to the length of the line associated with Q} 2 described above which for 
this example would be A S/2. Similarly the arc length associated with QJ 4 would be AS' U = R |+vt A0, ,/2. 
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Figure 3-3 Schematic of grid network and flow control area for discretization process 
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The area element dA will be made up of the parts of the flow control area In each of the four quadrants 
about the center point (l,j) and numbered based on the shaded half grid point that each contains. Thus 
AA 1 = (R( + Rj +V4 )A0 jA Sj/8, etc. and the discretized form of Equation (3-30) may be written as: 


AS 12 + Q 12 AS,2 + Q^AS^ + Q u AS; 4 - AS^ - Q** AS^ - Qa ASg - AS^, ■ ( 3- 31 ) 

+ P|)|(«ft + H r ) K i J ^AA 1 +(«6 +H r ),a J .iAA 4 + (aa ♦H l ) l j,^AA 3 + (a5 + H r ),_i.j,lAA 2 j| . 


The flow components on the left hand side of Equation (3-31) are obtained from discretization of Equations 
(3-21) and (3-22). A numbering system for the 9 pressures at the grid point (l,j) and the 8 surrounding 
points is shown in Figure 3-1, where P 1 * P, +1 j_ 1( P # ■ Py etc. The determination of the flows out of the 
sub-area containing the half grid point labeled 1 is discussed here as an example. The flow component 
Qjj, is determined from Equation (3-21). The derivative of the pressure normal to the line connecting points 
1 and 2 is evaluated at the intersection of the line with the grid line. The tangential derivative is evaluated 
at the half grid point, 1 as the average of the difference between P 3 and P 6 with that between P 2 and Pg, 
divided by AS,. Thus, 


lap 
r ae 


Pe-P» aP _ (Ps-PeWPa-Ps) 
R,A 0 , ’ as 2 AS, 


(for q; 2 ) . 


The flow component Q* 4 is determined in a similar manner from Equation (3-22). The normal derivative is 
approximated as the difference between P 2 and P 5 divided by AS, and the tangential derivative Is 
approximated as the average of the differences between P 3 and P z and P 6 and Pg, divided by R, +V4 A9j. 


3P 

as 


p a-P» lap 
as, ' r ae 


(P,-p 2 )*(P«-p«) 

2R.UA6, 


(for q; 4 ) . 


For both of the above flow components the pressure In the (1+P) term appearing In Equations (3-21) and 
(3-22) Is evaluated at the half grid point by averaging the four surrounding pressures (P 2 +P 3 +P 5 +Pg)/4. 
All of the remaining quantities (R, H r k v k* kg, k 4 , a , p and «) are evaluated directly at the half grid point 


The flow balances over the other quadrants are performed In a similar manner. For steady state conditions, 
the right hand side of Equation (3-31) will be 0 and the flow balance about any interior grid point (l,j) may 
be written in the form 
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F,(H r .P 1 .P 2 .P s .P 4 .P 5 .P # ,P 7 .P,.P 9 ) - 0 


(3-32) 


The definition of F.j may be extended to make Equation (3-32) applicable to the ends of the seal as well as 
the interior points by applying Equation (3-28) at the endpoints as follows: 


Fi, - P 5 - P, (1-1) and F mj - P, - P r (l-M) . 

The solution to Equation (3-32) may be used to provide all of the steady state quantities such as pressures, 
forces, moments, flow rate and power loss. The inclusion of the right hand side of Equation (3-31) will be 
necessary for determination of frequency dependent stiffness and damping coefficients which will be 
discussed later. 


Newton-Raphson linearization procedure 

The Newton-Raphson [5] procedure Is perhaps the most widely used method for obtaining solutions to non- 
linear systems of algebraic equations and is described in many textbooks on numerical methods such as 
Reference 5. A procedure similar to that used here is described In a paper by Artiles, Waiowit and Shapiro 
[6]- 

The procedure is started with an Initial pressure distribution that satisfies the end conditions given by 
Equation (3-28). A new set of approximations to the pressures in Equation (3-32), P k n * w may be obtained 
by linearizing Fg about a previously established set of approximations P k as follows: 



(3-33) 


where a forward difference 

3Fj ^ P|(Hr»Pl» — »P|t 4 ' t l»— »P>) ~ F|(H f ,Pf t -.,P 9 ) 

dP k T) 


may be used to numerically evaluate the partial derivatives. 

Pressures without the superscript new relate to the previous or "old'' approximation. It should be noted that 
the function Fg will not be 0 unless the pressures comprising its arguments are exact. If we go back to 
using grid notation for P (P,j in place of P 5 etc.) and introduce the column vector {Pj new } as the M new 
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pressures at the Jth column of grid points, Equation (3-33) may be written in the following form: 


| C J]{P,"**} ♦ |EJ]{P j T'} ♦ ID'HPjT) - (R 1 ) . 


(3-34) 


where 1C 1 ]. [E*] and [D>] are tri-diagonal matrices whose interior elements, from Equation (3-33), 


are 


Cu* ■ 4=*— . Ey-k “ , Dy** - ^ , k- -1,0,1 ; I -2,_,M-1 


i*fcj 




0P, 


1 * 0*1 


The interior elements of the column vector {R*} are 


“ E (C|!l*kP|*k,| + ^aPw^H * ^l!l*k^t*k.J*l) “ E|j . 

k *-1 

The above equations may also be applied to the comer elements to produce the result 

C i.i * Ci* - 1 . E* f1 - Eu, m - Di* t1 - Di.,, - 0 . R* - P, , Ri - P r . 

Equation (3-34) represents a linear system of simultaneous equations that may be solved by various matrix 
inversion procedures. The method used here Is the column or transfer matrix method, which is described 
in References 4 and 7. It has been used extensively in solving finite difference problems associated with 
various forms of the lubrication equations and produces accurate results In a fairly efficient manner. 
Convergence of the Newton-Raphson procedure is generally obtained within 3 - 6 iterations depending on 
degree of non-linearity and the accuracy required. 

Determination of loads, moments, torque and leakage 

The dimensionless loads and moments may be obtained by integrating the pressure distribution over the 
seal area as shown below: 


2* Sr 2 k S, 2, s, 

w x . / / Pcos0RdSd0, \fify * / / Psin0RdSd6, f j PRdSd0, 

OS, 0 8, OS, 


2k S, *, 8, 

M x m - j j Psin0RSdSd0, M y - / / Pco$0RSdSd0 . 
°s. os, 


(3-35) 
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The dimensionless torque Is obtained from integration of the shear stress given by Equation (3-26) over the 
seal area: 


2k Sr 

t-sign(u) J j xRdSdO . 
o s, 


(3-36) 


The sign(u) term has been added to make the torque positive when it opposes the net surface motion 
regardless of which surface (smooth or grooved) is moving. 


Finally, the dimensionless leakage flow, going into the seal at S * S, may be obtained from integration 
of C^, given by Equation (3-21), over the circumference of the seal: 


2 « 

Q h = / Q,Rd8 . (3* 37 > 

o 

The integrand in the above expression is evaluated by summing the flow components to the right of the first 
6 grid line in the same manner as that used in developing Equation (3-31). It should be noted that any 
value of S can be used since is Independent of S. 

The physical quantities corresponding to the dimensionless ones given above are 

M m - Rjfcfl*, . T.R||vpf .<««> 

In the above equations the loads W x and W y apply only to a cylindrical seal and W z applies only to a face 
seal. The leakage flow is the volumetric flow rate going Into the seal measured at pressure p ff 
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Determination of stiffness and damping coefficients 

Equation (3-20) with flow components and Q # given by Equations (3-21) and (3-22) represents a second 
order non-linear partial differential equation that may be used to define a second order non-linear operator 
G, such that 


G(P.H r ) - -4l(«»*H r )(1 + P)] . 
dt 


(3-39) 


The determination of P under steady state conditions, where the right hand side of Equation (3-39) is 0, was 
described earlier in this section. These steady state pressures will now be referred to as P. The various 
eccentricities and rotations used in determining H r from Equation (3-27) may, for convenience, be put in the 
form of a row matrix as 


[€*,♦.♦1. (face seal) 

4 

tl , (shaft seal) 


and Equation (3-27) may be written as 


H r * 1 + [€]{aj . 


(3-41) 


where the column vector {a}is given by 


(«} 


{-1 ,Ssin8,-Scos0} , (face seal) 

{-cos6, -sind ,S sine, -Scosd} , (shaft seal) 


(3-42) 


One could develop a perturbation analysis for prediction of stiffness and damping coefficients with the 
following procedure: (a) perturb say the ith component of the eccentricity matrix in Equation (3-40) by qe' , 
where q is a small parameter and c' Is time dependent; (b) express H r in the form H r = H + qe'{a> and 
the corresponding pressures as P = £ + q{P'}; (c) substitute the above expressions for P and H r in 
Equation (3-39); (d) expand the resulting expression neglecting terms of order q 2 and higher; (e) collect 
terms of order q . The resulting expression could be written in the form: 


Sf{P'} ♦ {bje' = -(«a - (1 + £){a}^ , 

dt dt 


(3-43) 


where Sf is a second order linear operator given by 
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(3-44) 


«■*■*•«* 


. A, A . i' 

^ae 2 




The coefficients in the above equations (Aj, Aj, (b) etc.) will depend on the coordinate variables as well 
as P, H and their various derivatives. Only the form of the above equations Is important to the numerical 
procedure under development and the significant amount of algebraic manipulation required to determine 
these coefficients will be shown to be unnecessary. 

If the time dependence of the eccentricity is restricted to oscillatory disturbances one may set c' = cP®* 
and look for solutions In the form {P'} ■ {P'je®®*, where {P*> Is complex but Independent of time. 
When this transformation is Introduced into Equation (3-43), the result Is 


$£{P*} ♦ (b) - -So[(«5 ♦ A)(P‘) ♦ (1 +P){a}] . t 3 - 45 ) 

The representation of {P*} and the eccentricity coefficients {a} as column vectors relates to the fact that 
each of the eccentricities must be perturbed to obtain the compete stiffness matrix but Equation (3-45) is 
solved Independently for each perturbation. The periodic boundary conditions given by Equation (3-29) also 
apply to Equation (3-39) (continuity of {P*} and 3{P*}/96 Is sufficient when H r and the spiral groove 
coefficients k 1r _,k 4 are continuous functions of 8 as they are here). The end boundary conditions given 
by Equation (3-28) become {P*} = 0 at S = S, and S = S r 

If Equation (3-45) were solved for {P*} subject to the above boundary conditions, all of the dimensionless 
stiffness and damping coefficients could be obtained by substituting {P } for P In Equation (3-35). The real 
parts of the computed forces and moments would be in phase with the eccentricity perturbations and 
constitute the dimensionless stiffness coefficients. Thus would correspond to the real part of W y 
computed from the component of {P*} associated with the perturbation In « x and S^ y would correspond 
to the real part of M x computed from the component of {P*} associated with the perturbation in e y etc. 
In a similar manner, the dimensionless damping coefficients which are 90* out of phase with the eccentricity 
perturbations would be obtained by dividing the Imaginary parts of the forces and moments computed In 
the manner described above by o. 

The parameter a , is a dimensionless disturbance frequency referred to as the ‘squeeze number* and is given 
by o * 2A.Q /« where Q is the angular velocity of the disturbance. The limiting form of the stiffness and 
damping coefficients as o - 0, is of interest as It applies to Incompressible flow, and the limiting stiffnesses 
are used in the homing procedure that has been Implemented for determining eccentricities from given loads 
which will be described later. This limiting form may be obtained by expressing Equation (3-45) in terms 
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of its real and imaginary parts as 


££{P») + (b) * o 2 («4+A){P # ) . 

(3-46) 

S£{P 8 ) --<«*+ A){P«}-(l+A){a> . 

(3-47) 


where 

{P*} - {P.} ♦ 3®{P,} . t 3 - 48 ) 

The column vectors {P,} and {P g } are the "stiffness* and "damping* pressures respectively. If one formally 
sets a = 0 In Equation (3-46) It decouples from Equation (3-47) and may be solved directly. Since the right 
hand side of Equation (3-46) becomes 0, the stiffness pressures are the same as those that would be 
obtained by computing the steady state pressures at a perturbed eccentricity, subtracting the unperturbed 
pressures and dividing by the eccentricity perturbation. This latter method is frequently used for computing 
steady state stiffnesses in incompressible flow and has been implemented here for the computation of 
"stiffnesses at 0 frequency" used in the above mentioned homing procedure. The 0 frequency damping 
pressures may be obtained by solving Equation (3-47) with {P g > as determined from the solution to 
Equation (3-46). 

The above discussion assumed that the perturbation coefficients In Equations (3-40) - (3-42) were 
determined prior to setting up the finite discretized equations for their solution. Identical results can be 
achieved by direct numerical perturbation of the difference equations. This approach, which has been 
Implemented here and is described below, avoids algebraic error in determining the perturbation coefficients 
and may be used in complex situations where analytical determination of the perturbation coefficients is not 
feasible. 

After desired convergence of the Newton-Raphson process has been achieved under steady (unperturbed) 
conditions one may denote the resulting steady state pressure vectors as {Pj} and the coefficient matrices 
as [ 5 1 ), etc. and Equation (3-34) may be written as 


+ IGWm) - {A 1 } - 


(3-49) 


One may now perturb the kth component of the eccentricity vector by an amount q , recalculate [fc*] at the 
new film thickness (but old pressure distribution, £) then subtract [6*J at the old film thickness and divide 
the difference by q to numerically obtain the derivative of (2jf| with respect to e k which will be denoted by 



[C lM ]. Thus 


tt u, ifi'ik., - 16 'iu 

The matrices [6* J ‘j, [D i,k J and {R**} are obtained in a simBar manner from the other coefficient matrices. 
If we introduce a disturbance to « k of magnitude e'ti , as was done in deriving Equation (3-43), then the 
change In the coefficient matrix [fc*] would bee'n [5 ,,k ] with corresponding changes In the other coefficient 
matrices. If we disturb Equation (3-49) by replacing {fy with + n {Pj k }. [£j] with [5*1 + e'n [5* ,k ], 
etc. and collect terms of order i \ , the following expression is obtained: 


If we set «' to unity In Equation (3-50) then <Pj k > will become the 0 frequency stiffness pressure (the 
change in steady state pressure per unit change In eccentricity). It should be noted that the coefficients of 
the 0 frequency stiffness pressures In Equation (3-50) are the same as those for the steady state pressures 
in Equation (3-49); only the right hand side has changed. Equation (3-50) thus represents the construction 
of the discretized form of Equation (3-43) when a ~ 0. In order to complete the process for a * 0, one 
may introduce the same disturbances to the right hand side of Equation (3-31), with H r « Fl + «'ii {a} and 
add the terms of order n to the right hand side of Equation (3-50). The terms to be added are 
-3([C J ]{P'j k }+{R*’ k }€')/8t, where [C*J are diagonal matrices whose components are 


C|-(a5 ♦ A)|»i j*a AA, ♦ (afi + ft) ( a i j_iAA 4 + (a3 ♦ + (a& ♦ft)|U,j*lAA 2 (3-51) 

and {R**} are column vectors whose ■ 


\ rnmonnontc : 


Hf* 


.(1 *P,)(a,^ | j4A 1 *a ) ^ | j4A 4 *a 1 \ ) jAA,-^ | j4A l )-(1 * . ( 3 . 52 ) 


The far right side of Equation (3-52) is a quadradlcally equivalent representation that was used in the 
computer program described In Section 3. One may now set e' * rP® 1 in Equation (3-50) and look for 
solutions in the form {Pj k } = {Pj k }®P*\ by introducing these substitutions into Equation (3-50) and 
combining terms to obtain the final set of linear difference equations for the complex stiffness pressures 


{P?}: 

IC-ihp; 11 ) M^HPh) ♦ ift'liPni -{R lk l -lC*l(^|l -t^nP'Hi • <3 ' 53) 
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where [C**] = [6*1 + SoIC 1 ] and {R i ’ k > = {R , ’ k > -3a{R i ’ k >. 


The system of equations given by Equation (3-53) has been solved by the column method In a directly 
analogous manner to that used in solving Equation (3-34). The principal difference lies in the fact that all 
of the matrix operations were performed using complex arithmetic. The dimensionless, frequency dependent 
stiffness and damping coefficients were computed from the complex stiffness pressures in the previously 
described manner. Relationships of the following type may be used to calculate the physical stiffness and 
damping coefficients from the dimensionless ones: 

K^-KoR*. K*-KoRo 2 R*. K^-KoRoR* t 3-54 * 

and 

« B 0 Ro 2 6 h . B*-B 0 Ro§* . (3_55) 

where Kq - p^/C and B 0 = 12pR 0 4 /C 3 . 

Optimization of groove parameters for maximum stagnation pressure in a concentric cylindrical seal 

Since spiral grooves are solely responsible for the axial stiffness of an aligned, gas lubricated face seal with 
parallel surfaces under steady state conditions, it is often desirable to optimize groove parameters for 
maximum axial stiffness. An optimization procedure for doing this has been implemented In the computer 
code SPIRALP described in Reference 8. The analogous situation is not as evident in a concentric gas 
lubricated cylindrical seal which will have considerable, if not maximum stiffness without spiral grooves. A 
large portion of the stiffness in the absence of spiral grooves wil be cross coupled, particularly at low values 
of A, thus giving rise to stability problems which may be alleviated with the use of spiral grooves. The 
criteria for optimizing groove geometry from a dynamic standpoint would thus depend on both the desired 
load capacity and the various other elements In the system affecting rotordynamlc performance. 

An alternate approach for developing a stand alone criterion for optimizing groove geometry in a cylindrical 
seal is to maximize the pressure gradient that the grooves can generate at stagnation. If the grooves are 
being used to pump against a pressure gradient, the maximum stagnation pressure gradient would represent 
the maximum pressure gradient that the grooves could pump against without allowing any net flow to go 
through. It would also represent the maximum axial pressure gradient that the grooves could generate in 
an aligned, symmetric herringbone bearing in the absence of an imposed pressure gradient In any event, 
the stagnation pressure gradient is a strong measure of spiral groove performance and even though 
optimizing It Is not a precise criterion for optimizing dynamic performance, computations obtained with 
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geometries optimised In this manner should provide a strong Indication of the maximum benefits obtainable 
with the use of spiral grooves. 


The stagnation pressure gradient for a cylindrical seal under concentric conditions may be obtained from 
Equation (3-22) by setting Q, = 0 (stagnation), dP/dQ = 0, H, = 1 (concentric). S = Z and R = 1 
(cylindrical seal). The resulting equation may be solved for A P/AZ making use of the definition of A 6 given 
by Equation (3-24) and the definitions of k, and k 4 given by Equation (3-25) to obtain the following 

relationship 


JP „ AS ~«> s| npco3P(I«-11 
aZ «(1 -«)(!* -Ifsln 2 ? +1* 

The right hand side of the above equation may be treated as a function of a , P and i F = 1 + « ) and has 
a maximum value of dP/dZ = 0.09118 AS at = 0.5, P^ » 0.2736 (15.68? ) and 5 « 2.653. 

The variation of the pressure gradient near the optimum point is shown In Figure 3- 4. The curve marked 
a was obtained by holding p and 5 at their optimum values and varying «. The other curves were obtained 
in an analogous manner. The curves show the sensitivity of the optimum pressure gradient to the various 
parameters and verify the existence of a relative maximum at the optimum point 

Other approaches to the optimization problem are given In Reference 2 for spiral groove bearings and 
Reference 9 for spiral groove viscous pumps. 
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3.2 Description ot Computer Code SPIRALS and Subroutine SPIRAL 

A FORTRAN subroutine. SPIRAL has been written to Implement the analysis developed under Section 3.1.2 
In dimensionless form. The analysis has been programmed In this tarn, to pennlt easy Incorporation Into 
me knowledge base system currently under development SPIRAL and «s associated sub-programs 
constitute a self contained system that has no Input-output other ttan the arguments passed to It through 
SPIRAL and Is thus Independent of the operating system. SPIRAL has bear compiled, m Its present torn, 
with Version 5.0 of the MlcrosofP Fortran Compiler and should work with many other compilers with 
relatively little modification. Significant user Information Is Included In the Users Manual. Reference (10). 

The analytical procedure contained In Section 3-2 has been oriented toward determining pressure 
distribution, load. How. torque, stiffness and damping for a given Urn thickness distribution. In practice It 
Is often desirable to determine the equilibrium film thickness or eccentricities from prescribed loads and 
possibly moments. SPIRAL provides a homing option for determining the eccentricities based on the steady 
state bearing stiffnesses. This homing option Is based on the procedure described below. 

If one were to write the dimensionless load and eccentricity as column vectors (W)and(r) (transpose row 
matrix M) and take the previous estimate (or Initial guess) of« as M* and the load vector computed 
from (cloy as (W),*,, the steady state stiffness matrix [fcl cor*l be used to arrive at a new approximation 
for «.) The method for doing this Is shown by first writing the equation for the change In load as 
<W> - - M<(«>- («W The new approximation to *} Is obtained by Inverting the stiffness 

matrix and solving for (r> as (<> • * (Kl’tW - {*W This a P proach 13 ^ ^ "* 

application of the Newton-Raphson method for determining the eccentricities. 


Whie the above approach can be very effective It can also diverge » the Initial guesses are bad. This 
divergence Is usually accompanied by the generation of negative film thicknesses In the course of the 
Iteration process. In order to aaempt to coned this problem, an optional numerical damping algorithm has 
been Implemented which replace, M with *> - M* ♦ - <*W when th. originally 

calculated value of {e} would result in a negative fflm thickness. 

The cell method of discretization is designed to obtain quadratic accuracy. Numerical testing Indicates that 
this has apparently been achieved. One may make use of this property to obtain greater accuracy, (or the 
same degree of accuracy with coarser grids and ensuing reductions In computer time) with the use of 
Romberg extrapolation. Suppose for example we computed the dimensionless torque T with a coarse grid 
and denoted It by t c then halved the grid spacing In both directions and recomputed t denoting it as % 
(subscript denotes fine grid). If the truncation error were to approach 0 as the square of the grid spacing 
and T r were the true solution then <i f -T r ) * <t c - 1,)/4. or t, = (4T f - tj/3. The above extrapolation 
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can. in principal, increase the rate of convergence from quadratic to cubic. The subroutine SPIRAL, provides 
the option of implementing Romberg extrapolation. 

The logic used in SPIRAL for performing the pressure Iterations, computing stiffness and damping 
coefficients, homing in on eccentricities and implementing Romberg extrapolation is shown in Figure 3-5. 
It can be seen there that when the homing process Is Implemented. It is completed for both coarse grid 
and fine grid solutions prior to performing the Romberg extrapolation. The extrapolation is thus performed 
with solutions obtained at two different displacements. When the displacements are specified, extrapolations 
are performed with solutions obtained at the same displacement, which is believed to be a more accurate 
approach. If one were to compute displacements for a given loading and then recompute the loading from 
the displacements using Romberg extrapolation for both computations the computed loading would thus 
differ slightly from the input loading even though all tolerances were met The degree of this difference will 
depend on the grid size and caution should be exercised In using Romberg extrapolation when homing in 
on the displacements with very coarse grids. 
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Figure 3-5 Flow diagram for logic used in SUBROUTINE SPIRAL 














3.3 Sample Problems 


A number of sample problems have been prepared to demonstrate the behavior and various features ol the 
computer program. They are Intended primarily tor lustration and do not necessarily represent 

recommended seal designs. 

Oases t - 3 serve to show the Improvement In accuracy that can be obtained with the use of Romberg 
extrapolation tor a concentric, asymmetric cylindrical seal The a* I. divided Into two regions ol equal 
length as shewn In Figure 3-7. The stationary surface m the Urst legion has a groove geometry optimized 
lo, maximum stagnation pressure as desedbed at the end ol Section 3-2. The grooves are oriented to 
produce a pumping component In the positive axial direction to partially offset the larger negative one 
caused by the Imposed pressure gradient The second region Is smooth. The trase t results were obtained 
without the use of Romberg extrapolation. Romberg extrapolation was used In Case 2 with coarse grid 
solution obtained lor the same grid geometry as tisad In Case 1. Since the grid spring Is halved m each 
direction when obtaining the line grid solution. Case 2 should represent a much more accurate solution than 
Caset. It also took approximately 1 1 times as long to run. Romberg extrapolation was used In Case 3 with 
twice the grid spacing as that used m Case 2 and took ody 25* longer to run than Case 1. Thedlrect 
stiffness coefficients, K„, calculated lo, Cases 1 • 3 are 54684. 57622 and 57342 Ib/m, respectively. Using 
Case 2 as a standard, the error In the Case 1 stillness Is 5 . 1 * while the error In the Case 3 stiffness Is only 

.5%. 

A symmetric ‘herringbone groove’ pattern Is used In Case 4 with the same overall geometry as that used 
in Cases 1 - 3. The groove pattern on the stator in region 2 Is the same as that for region 1 with the 
exception of the sign of the groove angle. The grid geometry Is the same as that used in Case 1 but the 
operating conditions differ in that there Is no Imposed pressure gradient and the shaft Is displaced In the 
x direction and tilted about the y axis. It can be seen that the Imposed displacement and tDt produces non- 
zero values for the calculated forces and moments. The results of Case 5 were obtained by prescribing the 
forces and moments computed for Case 4. The initial guess for the shaft displacement for Case 5 was taken 
to be somewhat larger than prescribed for Case 4 and Initial guess for the tnt was taken as 0. The 
displacements calculated for Case 5 are essentially the same as those Imposed In Case 4. Case 6 Illustrates 

the use of the program with SI units. 

The remaining 2 cases correspond to a mechanical face seal under a very high imposed pressure gradient 
with spiral grooves on the outside surface of the stator oriented to pump inward with the pressure gradient 
The stator geometry is shown schematically in Figure 3-8. with the land width somewhat enlarged for clarity. 
The inward pumping is Induced by the counterclockwise motion of the rotor. Solutions to this problem will 
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be approximate in nature In that choking Is likely to occur (not treated here as a result of assumed 
Isothermal How with negligible Inertia) which w» raise the eftectlve Bm pressure at the Inside radrus to a 
value somewhat higher than mat prescribed. These cases are provided to llusbate the use of the program 
with a lace seal and the evaluation ot the Internal accuracy ol the program. Case 8 was obtained by halving 
the grid spacing used In Case 7, In both directions. This procedure provides a test for truncation error, 
which Is small m this case, that Is recommended lor frequent use m determining appropriate grid spacing. 
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Figure 3-6 Schematic of shaft seal for Cases 1-3 
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pumping against pres. grad. 


( CASE 1 ) Asymmetric cyl. seal with grooves 
SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 5.0000E-04 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.0000E+04, 1.0000E+04 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.4700E+01. 6.4700E+01 (PSI) 

VISCOSITY- 2 . 9000E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 3 0 

CALCULATED FORCK IN X.Y DIRECTIONS - gRo2E-15 2927E- 14* 1 (IN-Lb! 
CALCULATED MOMENTS ABOUT X,Y AXES - -8.600ZE-.L3, *.***/*• * \ 

MINIMUM FILM THICKNESS - 5.0000E-04 (IN) 

FLOW - -2. 1757E+00 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE- 7.4543E-02 (IN-LB), FILM POWER LOSS - 1.1827E-02 (HP) 

COMPRESSIBILITY NUMBER - 4.9582E+00, SQUEEZE NUMBER - 9.9163E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphi 

Bpsi 


5 46&eIo 4 3.?9®04 -IfUfE&A 

VoQlfiE+04 5 4684E+04 6.0468E+04 -2.7046E+04 

-8 1693E+03 -3 4030E+03 1.0593E+04 

8 16931+03 7 8082E+03 -l!0593E+04 -3.4030E+03 

"l 1158E+02 -6:2395itoi HlO^. 

•fiiH i:«? 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB -SEC 

IN-LB-SEC 

IN-LB-SEC 
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( CASE 2 ) Romberg extrapolation with coarse grid the same as Case 1 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 5.0000E-04 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.0000E+04, 1.0000E+04 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.4700E+01, 6.4700E+01 (PSI) 

VISCOSITY - 2.9000E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 3 0 

CALCULATED FORCES IN X,Y DIRECTIONS - 3.0897E-13. -^7781E-13 (13) 

CALCULATED MOMENTS ABOUt X,Y AXES - -9.7768E-14, 1.3279E-13 (IN LB) 

MINIMUM FILM THICKNESS - 5.0000E-04 (IN) 

FLOW - -2.1749E+00 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE - 7.4550E-02 (IN-LB), FILM POWER LOSS — 1.1829E-02 (HP) 

COMPRESSIBILITY NUMBER - 4.9582E+00, SQUEEZE NUMBER - 9.9163E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kx 

Kphi 

& Si 

By 

Bphi 

Bpsi 


5 ySio* 3 JsWBo* 

-3 2577E+04 5.7622E+04 6 4855E+04 -2.9332E+04 

8 1664E+03 -8.4866E+03 -3.6583E+03 1.2003E+04 

8 - 4866E+03 8.1664E+03 -1.2003E+04 -3.6583E+03 

i * 1 -5 1897E+01 z 0183E+01 6.7544E+01 

\ 9 1897E+01 1 * 1837E+02 -6 . 7 544E+01 2 . 0183E+01 

-1 2229E+01 -1 ! 8217E+01 1.9551E+01 

l!8217E+01 -1 . 2229E+01 9.0828E+00 1.9551E+01 


FORCE UNIT 

LB 

LB 

IN-LB 
IN-LB 
LB- SEC 
LB- SEC 
IN-LB-SEC 
IN-LB-SEC 
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( CASE 3 ) Romberg extrapolation with fine grid the same as Case 1 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 5.0000E-04 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.0000E+04, 1.0000E+04 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.4700E+01. 6.4700E+01 (PSI) 

VISCOSITY- 2.9000E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 3 0 

MINIMUM FILM THICKNESS - 5.0000E-04 (IN) 

FLOW - -2.1750E+00 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE- 7.4549E-02 (IN-LB), FILM POWER LOSS - 1.1828E-02 (HP) 

COMPRESSIBILITY NUMBER - 4.9582E+00, SQUEEZE NUMBER - 9.9163E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphi 

Bps! 


5 ?3®04 3.L9eI<) 4 -5 h WlE&, 

-32559E+04 5.7342E+04 6.4238E+04 -2.9001E+04 

812990E+03 -8!5936E+03 -3.7301E+03 
ft ^93fiF+03 8 2990E+03 -1.1766E+04 -3.7301E+03 

l ' 1 7fifiF+02 -3 1507E+01 1 9728E+01 6.6992E+01 

3 * 15071+01 i:i766E+02 -6:6992E+01 1.9728E+01 

11432E+01 -1 8433E+01 1.9494E+01 -8.7410E+00 

l! 8433E+01 .i:i932ESl 8:7410E+00 1.9494E+01 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB-SEC 

IN-LB-SEC 

IN-LB-SEC 
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( CASE 4 ) Displaced and tilted symmetric cyl. seal, no pres. grad. 

SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 5.0000E-04 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.0000E+04, 1.0000E+04 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.4700E+01, 1.4700E+01 (PSI) 

VISCOSITY- 2 . 9000E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 4 0 

spawn rn&W&SBF 

MINIMUM FILM THICKNESS - 2.5000E-04 (IN) 

FLOW - 2.0302E-02 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE- 5.8726E-02 (IN-LB) , FILM POWER LOSS - 9.3178E-03 (HP) 


COMPRESSIBILITY NUMBER - 4.9582E+00, SQUEEZE NUMBER - 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphi 

Bps! 


* fcftliM -f S 85^2 

SiilaSiiol i:f8gK8j 

i 7fi4oF+02 7 3977E+01 -1.4889E+03 4 . 8095E+03 

*9 " 3280E+01 -l‘ 3471E+01 1 0865E-01 6.3426E-01 

I’^lnortni 22299E+01 -1 3412E-01 2.1028E-01 

7 * 7952E^02 -1 0697E-01 4i5324E+00 -1.0203E+00 

4.7944E-01 -7. 1164E-02 9.5944E-01 4.8461E+00 


9 . 9163E+00 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB -SEC 

LB-SEC 

IN-LB-SEC 

IN-LB-SEC 
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( CASE 5 ) Same seal as Case 4 with applied forces Instead of dlspl. 
SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00. 5.0000E-04 ( ) 

ROTATION SPEED. DISTURBANCE SPEED - 1.0000E+04, 1.0000E+04 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.4700E+01, 4700E+01 (PSI) 

VISCOSITY- 2.9000E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ERROR CODE - 0, ITERATIONS IN HOMING PROCESS - 3 

SMI8 ?RTS A AB^iFAX£S Y - DI 9 E 8499E?10, ' <*»> ' 

ITERATIONS IN LAST PRESSURE CALCULATION - 1 

gftgg&B SiMs%oik Y xT^s°!2 s 2.62t?l?Sf!°°J.U«Sl-8! t0 ?iN®^ 

MINIMUM FILM THICKNESS - 2.5000E-04 (IN) 

FLOW - 2 . 0301E-02 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE- 5 . 8726E-02 (IN-LB), FILM POWER LOSS - 9.3178E-03 (HP) 

COMPRESSIBILITY NUMBER - 4.9582E+00, SQUEEZE NUMBER - 9.9163E+ 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 

5IS p - , HAVSIaa Q Sa® 03 J^SE^l 2*“ 


•lii mi mi f is 

Hi m Hi am 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB -SEC 

LB-SEC 

IN-LB-SEC 

IN-LB-SEC 
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( CASE 6 ) Same as Case 4 with SI units 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 5.0800E-02, 5.0800E-02, 1.2700E-05 (m) 

ROTATION SPEED, DISTURBANCE SPEED - 1.0000E+04, 1.0000E+04 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.0135E+05, 1.0135E+05 (Pa) 

VISCOSITY - 1.9995E-05 (Pa-SEC), AMBIENT PRESSURE - 1.0135E+05 (Pa) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 4 0 

shm Ss« Y xTsi^ s 

MINIMUM FILM THICKNESS - 6.3500E-06 (m) 

FLOW - 3.3268E-07 (m**3/SEC) MEASURED AT 1.0135E+05 (Pa) 

TORQUE- 6.6352E-03 (N-m) , FILM POWER LOSS - 6.9483E+00 (WATT) 

COMPRESSIBILITY NUMBER - 4.9584E+00, SQUEEZE NUMBER - 9.9167E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


7 feiSW iJoftW -! h 70i0E?i2 

-1 6523E+06 7.0953E+06 1.1274E+03 5 . 7639E-i^2 

. 1 * 3110E+02 9 1 5888E+02 5 . 4708E+02 1 . 7032E+02 

-5 6446E+02 3.2904E+02 -1.6822E+Q2 5.4341E+02 

Zt # n7fifiF+03 -2 3592E+03 4*8331E-01 2.8213E+00 

7 * 3833E+03 3 9051E+03 -5.9660E-01 9.3536E-01 

-3 2450E-01 -4.7582E-01 5.1209E-01 -1.1528E-01 

2!l326E+00 -3.1658E-01 1.0841E-01 5.4754E-01 


FORCE UNIT 

N 

N 

N-m 

N-m 

N-SEC 

N-SEC 

N-m-SEC 

N-m-SEC 
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Figure 3-7 Stator with inward pumping grooves for Cases 7 and 8 
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( CASE 7 ) Face seal for pipe line compressor 


SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS - 3.7930E+00, 4.5050E+00, 1.0000E-04 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.4500E+04, 1.4500E+04 (RPM) 

INSIDE, OUTSIDE PRESSURE - 1.4700E+01, 9.1470E+02 (PSI) 

VISCOSITY - 1.7500E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 4 0 

CALCULATED FORCE IN Z DIRECTION - 3j 4246E+03 (LB) , TW T n \ 

CALCULATED MOMENTS ABOUT X,Y AXES - -1.0069E-12, -5.7018E-12 (IN-LB) 

MINIMUM FILM THICKNESS - 1.0000E-04 (IN) 

FLOW - -1.3026E+02 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE- 3 . 2695E-01 (IN-LB), FILM POWER LOSS - 7.5219E-02 (HP) 

COMPRESSIBILITY NUMBER - 5.5030E+02, SQUEEZE NUMBER - 1.1006E+03 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kz 

Kphi 

Kpsi 

Bz 

Bphi 

Bps! 


7.SlfflPo« -ftMSBio -? S WJ% 01,11 

1-881:8 t:88Kg flfcft, 

9 1 7180E-04 1 . 7161E-04 LB-SEC 

-2i2269E-07 6.6669E+02 -7.0380E+01 

-1.0207E-06 7.0380E+01 6.6669E+02 IN-LB-SEC 
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( CASE 8 ) Same as case 7 with half the grid spacing in each direction 
SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD. REFERENCE FILM THICKNESS - 3.7930E+00, 4.5050E+00, 1.0000E-04 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.4500E+04, 1.4500E+04 (RPM) 

INSIDE, OUTSIDE PRESSURE - 1.4700E+01, 9.1470E+02 (PSI) 

VISCOSITY- 1.7500E-09 (PSI-SEC), AMBIENT PRESSURE - 1.4700E+01 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 4 0 

3 5 4 §f 2 2!E°j2. (L ! ) 5739E-12 (IN-LB) 

MINIMUM FILM THICKNESS - 1.0000E-04 (IN) 

FLOW - -1.3026E+02 (IN**3/SEC) MEASURED AT 1.4700E+01 (PSI) 

TORQUE- 3.2695E-01 (IN-LB), FILM POWER LOSS - 7.5219E-02 (HP) 

COMPRESSIBILITY NUMBER - 5.5030E+02, SQUEEZE NUMBER - 1.1006E+03 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 
nT qp 2 (IN) phi (RAD) psi (RAD) FORCE UNIT 

g- :rt ssi : li 

fe[ iiSSH l:8M Wm 
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3.4 Verification 


SPIRALG has been compared with the results of two computer codes. The first of these Is MTI Computer 
Code PN471 which Is fully described in Reference 3. The program is based on perturbation analyses and 
is only applicable to bearings and seals operating In the concentric, aligned position fe x *« y “4> =♦=<>). 
The program does not predict loads and moments that would occur at finite displacements but It does 
predict stiffness and damping values as well as flow, torque, and power loss for spiral groove bearings as 
well as cylindrical and face seals. Since the program solves ordinary rather than partial differential equations 
It can be made to rapidly produce highly accurate results for evaluating the accuracy of SPIRALG. The 
second MTI computer code, named GASBEAR. Is used to verify SPIRALG under displaced and misaligned 
conditions. GASBEAR was written for use in conjunction with plane journal bearings and cylindrical seals. 

It does not treat spiral grooves or face seals. Since SPIRALG does not contain any special Instructions for 
treating concentric behavior and has relatively few instructions for distinguishing between face and cylindrical 
seals, the above two programs should provide reasonable verification. Since the treatment of the effects 
of spiral grooves under eccentric and misaligned conditions Is believed to be new. terms that become 
significant only under those conditions remain unverified. 

The results of 7 verification tests are reported on the following pages. A SPIRALG output listing followed 
by the relevant output from the verification code, converted to equivalent units and format. Is given for each 
case. The somewhat strange looking input values (unit ambient pressure, high RPM but low viscosity etc.) 
were selected to simplify the conversion process between dimensional quantities and the dimensionless 
ones that were used throughout the development of the code. The compressibility number of A « 10 used 
for all cases and the seal pressure ratio of 2 used for imposed pressure gradients should be typical of many 
practical applications under fairly compressible conditions. 

Cases 1 - 6 show comparisons between SPIRALG and PN471. Romberg extrapolation was used for each 
of these cases, with 21 grid points In the circumferential direction and 5 sub-intervals In each of the two 
regions in the transverse direction for the coarse grid solution. The fine grid solutions thus use 41 
circumferential points and 10 sub-intervals per region in the transverse direction. 

The first case verifies stiffness and damping values for a synchronous disturbance acting on a cylindrical 
seal with a herringbone groove pattern and an Imposed pressure gradient. Cases 2 and 3 verify the 
differences in stiffness and damping values predicted to occur for a cylindrical seal when grooves are placed 
on the rotor (with groove angles reversed) rather than on the stator. The static quantities (flow, torque and 
power loss) remain unchanged. Cases 4 - 6 show comparisons between SPIRALG and PN471 for a spiral 
groove face seal with an imposed pressure gradient at three different disturbance frequencies; zero (case 
4), synchronous (case 5) and ten times synchronous (case 6). 
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Case 7 shows a comparison between SPIRALG and GASBEAR for and eccentric, tilted cylindrical seal with 
an Imposed pressure gradient. The grid size was chosen to match the maximum size allowable for the 
available version of GASBEAR. A separate program was written to perform Romberg extrapolations with 
the results of GASBEAR. The agreement between the two programs is good. The apparent discrepancy 
between the moments about the y axis is a result of the fact that the component Is very small. The relative 
error obtained by dividing the discrepancy by the absolute magnitude of the moment vector Is 0.33%. 
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( CASE 1 ) Concentric cyl. seal with pres. grad, lamda-10, sigma-20 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, 1.9099E+05 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 2.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY- 8 . 3333E-11 (PSI-SEC), AMBIENT PRESSURE - 1.0000E+00 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 2 0 

CALCULATED FORCES IN X,Y DIRECTIONS - ;li0902E-15 . ;®aU 5 ?|" 1 § tn Q‘S} 
CALCULATED MOMENTS ABOUf X,Y AXES - -9.3565E-16, -2.4405E-16 (IN-LB) 

MINIMUM FILM THICKNESS - 1.0000E-03 (IN) 

FLOW - 1.2816E+01 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

TORQUE- 1. 7015E-02 (IN-LB), FILM POWER LOSS - 5.1562E-02 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - 2.0000E+01 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kx 

Kphi 

Kpsi 

By 

Bphi 

Bpsi 


5 . u $&> 3 Yhimk 

:?:81”lt8i -imt®. jB8 ' : f 

1 7773E+01 -1.0824E+02 - 1 . 9521E+02 7.0561E+02 

7 ! 3200E-02 -5. 5832E-02 -1.5172E-02 -2.9538E-02 
5 . 5832E-02 7.3200E-02 2.9538E-02 -1-5172E-02 

1.9442E-03 6 . 9321E-03 2.7091E-02 -1.1690E-02 

-6.9321E-03 1 . 9442E-03 1.1690E-02 2.7091E-02 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB-SEC 

IN-LB-SEC 

IN-LB-SEC 
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COMPARISON OF CASE 1 WITH PN471 


FLOW - 
TORQUE 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphi 

Bps! 


1 . 282E+01 (IN**3/SEC) MEASURED AT 1.0000E+00 (PSI) 

- 1 . 701E-02 (IN-LB), FILM POWER LOSS - 5.156E-02 (HP) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


5.46&8E+03 
-9.0292E+01 
-1.0809E+02 
1 . 7525E+01 
7.3215E-02 
5 . 5805E-02 
1.9403E-03 
-6.9195E-03 


9 .$2^2E+01 
5.4608E+03 
-1.7525E+01 
-1.0809E+02 
-5.5805E-02 
7 . 3215E-02 
6.9195E-03 
1.9403E-03 


phi (RAD) 
-E . 8846E+01 
-2.6823E+02 
7.0356E+02 
-1.9537E+02 
-1.5192E-02 
2 . 9534E-02 
2.7097E-02 
1. 1667E-02 


psi (RAD) 
2.6823E+02 
-6.8846E+01 
1.9537E+02 
7.0356E+02 
-2.9534E-02 
-1.5192E-02 
-1.1667E-02 
2.7097E-02 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB- SEC 

IN-LB-SEC 

IN-LB-SEC 
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( CASE 2 ) Concentric asymmetric cyl. seal, lamda-10, sigma-0 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - l.OOOOE+OO, 1.0000E+00 (PSI) 

VISCOSITY- 8 . 3333E-11 (PSI-SEC), AMBIENT PRESSURE - l.OOOOE+OO (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 2 0 

CALCULATED FORCES IN X.Y DIRECTIONS - -2.8151E-15, -1.6233E-15 (LB) 
CALCULATED MOMENTS ABOUt X,Y AXES - 2.7524E-16, -5.2158E-17 (IN-LB) 

MINIMUM FILM THICKNESS - 1.0000E-03 (IN) 

FLOW - 3 . 5000E+00 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

TORQUE- 1.8850E-02 (IN-LB), FILM POWER LOSS - 5.7122E-02 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - 0.0000E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 
x 


DISP. 

Kx 

Kphl 

Sr 1 

By 

Bphi 

Bpsi 


y (IN) phi (RAD) 
1.5689E+03 4.996QE+62 


4. 96^7E+03 

-1.5689E+03 4.9657E+03 -2.6673E+02 
5.3012E+01 -3 . 6877E+02 6.0778E+02 
3. 6877E+02 5.3012E+01 -5.3984E+02 
-6.9953E-02 -1.1123E-01 4.3933E-02 
1. 1123E-01 -6 . 9953E-02 1.2462E-02 
-8. 3037E-03 -4.4522E-03 2.2471E-02 
4.4522E-03 -8.3037E-03 4.3624E-02 


psi (RAD) 

2 . 66/3E+02 
-4.9960E+02 
5 . 3984E+02 
6.0778E+02 
-1.2462E-02 
4.3933E-02 
-4. 3624E-02 
2.2471E-02 


FORCE UNIT 

LB 

LB 

IN-LB 
IN-LB 
LB -SEC 
LB -SEC 
IN-LB-SEC 
IN-LB-SEC 
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COMPARISON OF CASE 2 WITH PN471 


FLOW - 
TORQUE 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphi 

Bps! 


3 . 500E+00 (IN**3/SEC) MEASURED AT 1.0000E+00 (PSI) 

- 1 . 885E-02 (IN-LB), FILM POWER LOSS - 5.712E-02 (HP) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


x (IN) y (IN) 

4.9648E+03 1.5692E+03 

-1 . 5692E+03 4.9648E+03 

5.3156E+01 -3.6793E+02 
3 . 6793E+02 5.3156E+01 

-6 . 9955E-02 -1.1143E-01 
1. 1143E-01 -6 . 9955E-02 
-8.3373E-03 -4.4940E-03 
4.4940E-03 -8.3373E-03 


phi (RAD) 
-4. 9982E+02 
-2.6552E+02 
6 . 0602E+02 
-5.3900E+02 
4. 3975E-02 
1.2344E-02 
2 . 2617E-02 
4. 3536E-02 


psi (RAD) 
2.6552E+02 
-4.9982E+02 
5 . 3900E+02 
6.0602E+02 
-1. 2344E-02 
4.3975E-02 
-4.3536E-02 
2.2617E-02 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB -SEC 

IN-LB-SEC 

IN-LB-SEC 
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( CASE 3 ) Same case with grooves on moving surf. f groove angle rev. 
SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS GROOVED 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, 0.0000E+00 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 1.0000E+00, 1.0000E+00 (PS I) 

VISCOSITY - 8.3333E-11 (PSI-SEC), AMBIENT PRESSURE - 1.0000E+00 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 2 0 

CALCULATED FORCES IN X.Y DIRECTIONS - ;342 2 ?f' 16 i 9 ^sIf 5 ? f^riN^i 
CALCULATED MOMENTS ABOUT X,Y AXES - 1.6213E-16, 6.238-6E-17 (IN-LB) 

MINIMUM FILM THICKNESS - 1.0000E-03 (IN) 

FLOW - 3 . 5000E+00 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

TORQUE- 1.8850E-02 (IN-LB), FILM POWER LOSS - 5.7122E-02 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - 0.0000E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


4.52^1E+03 1.^9^3E+03 -S^TE+jta 5?jta£|js+i)2 

•UlUlfol R81 1:t8i|R8§ 1: H 

4 7068E+02 3.4694E+02 -5.1362E+02 5.9016E+02 

-2.3432E-02 -1.2667E-01 3.4823E-02 -§-?373E-03 

l!2667E-01 -2.3432E-02 5.9373E-03 3.4823E-02 

-2 . 7769E-02 6.0874E-03 2.3559E-02 -^-^038E-02 

-6.0874E-03 -2.7769E-02 4.4038E-02 2.3559E-02 


FORCE UNIT 

LB 

LB 

IN-LB 
IN-LB 
LB -SEC 
LB -SEC 
IN-LB-SEC 
IN- LB -SEC 
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COMPARISON OF CASE 3 WITH PN471 


FLOW - 
TORQUE 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphl 

Bps! 


3 . 500E+00 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

- 1 . 885E-02 (IN-LB), FILM POWER LOSS - 5.712E-02 (HP) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


4. 52^6E+03 
-1.2934E+03 


1.^9^4E+03 
4. 5236E+03 


3.4697E+02 -4.6939E+02 
4. 6939E+02 3.4697E+02 
-2.3471E-02 -1. 2673E-01 
1.2673E-01 -2 . 3471E-02 
-2.7798E-02 5.9155E-03 
-5.9155E-03 -2.7798E-02 


phi (RAD) 
-3.6992E+02 
-2.1883E+02 
5 . 8843E+02 
-5. 1292E+02 
3.4849E-02 
5. 8520E-03 
2.3693E-02 
4. 3962E-02 


psi (RAD) 
2.1883E+02 
-3.6992E+02 
5 . 1292E+02 
5 . 8843E+02 
-5.8520E-03 
3.4849E-02 
-4.3962E-02 
2.3693E-02 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB-SEC 

IN-LB-SEC 

IN-LB-SEC 
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< CASE 4 ) Face seal, no tilt, with pres, grad., lamda-10, sigma-0 
SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS - 1.0000E+00, 2 . 0000E+00 , 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, 0.0000E+00 (RPM) 

INSIDE, OUTSIDE PRESSURE - 2.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY- 8 , 3333E-11 (PSI-SEC), AMBIENT PRESSURE - 1.0000E+00 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 2 0 

CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES - -3.0782E-16, 1.1794E-15 (IN-LB) 

MINIMUM FIIM THICKNESS - 1.0000E-03 (IN) 

FLOW - 2.2763E+01 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

TORQUE - 2.2645E-03 (IN-LB), FILM POWER LOSS - 6.8621E-03 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - O.OOOOE+OO 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


„ z (?N) phi (RAD) psi (RAD) FORCE UNIT 
3 . 9524E+02 -4.8585E-04 -4.8576E-64 LB 

-6 . 1349E-07 2.0998E+02 8.2256E+01 IN-LB 

-3.0546E-07 -8.2256E+01 2.0998E+02 IN-LB 

2 . 9491E-02 -7 . 3530E-07 -7.3551E-07 LB-SEC 

-1. 1280E-09 7.0914E-03 -1.8082E-03 IN-LB-SEC 

-2. 3148E-09 1 . 8082E-03 7.0914E-03 IN-LB-SEC 
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COMPARISON OF CASE 4 WITH PN471 


CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 

FLOW - 2 . 2763E+01 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

TORQUE - 2 . 2644E-03 (IN-LB), FILM POWER LOSS - 6.8620E-03 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kz 

z (IN) 

3. 9524E+02 

phi (RAD) 

psi (RAD) 

FORCE UNIT 
LB 

Kphi 


2.0998E+02 

8 . 2240E+01 

IN-LB 

Kpsi 

2.9485E-02 

-8 . 2240E+01 

2.0998E+02 

IN-LB 

Bz 



LB-SEC 

Bphi 


7.0901E-03 

-1.8058E-03 

IN-LB- SEC 

Bps! 


1. 8058E-03 

7.0901E-03 

IN-LB-SEC 


(HP) 
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( CASE 5 ) Face seal, no tilt, with pres, grad., lamda-10, sigma-20 
SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS - l.OOOOE+OO, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, 1.9099E+05 (RPM) 

INSIDE, OUTSIDE PRESSURE - 2.0000E+00, l.OOOOE+OO (PSI) 

VISCOSITY- 8 . 3333E-11 (PSI-SEC), AMBIENT PRESSURE - 1.0000E+00 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 2 0 

CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES - -3.0782E-16, 1.1794E-15 (IN- LB) 

MINIMUM FILM THICKNESS - 1.0000E-03 (IN) 


FLOW - 

2.2763E+01 

(IN**3/SEC) 

MEASURED AT 

l.OOOOE+OO 

(PSI) 

TORQUE 

- 2.2645E-03 

(IN-LB), 

FILM POWER 

LOSS - 6.8621E-03 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - 

2.0000E+01 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 


DISP. 

. z (IN) 

phi (RAD) 

psi (RAD) 

FORCE UNIT 


Kz 

5 . 3566E+02 

-2.9756E-04 

-2 . 9746E-04 

LB 


Kphi 

-2 . 1294E-07 

2 . 4013E+02 

7 . 0968E+01 

IN-LB 


Kpsi 

2 . 5103E-07 

-7 . 0968E+01 

2.4013E+02 

IN-LB 


Bz 

2 . 7697E-02 

3 . 8068E-09 

3 . 8070E-09 

LB-SEC 


Bphi 

2 . 2173E-12 

6 . 7533E-03 

-1 . 6346E-03 

IN-LB-SEC 


Bpsi 

-2.0783E-12 

1 . 6346E-03 

6. 7533E-03 

IN-LB-SEC 
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COMPARISON OF CASE 5 WITH PN471 


CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 


FLOW - 

2.2763E+01 

(IN**3/SEC) 

MEASURED AT 

1 . OOOOE+OO (PSI) 

TORQUE 

- 2 . 2644E-03 

(IN-LB), 

FILM POWER 

LOSS - 6.8620E-03 (HP) 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / DISP. UNIT ) 

DISP. 

Kz 

z (IN) 

5 . 3557E+02 

phi (RAD) 

psi (RAD) 

FORCE UNIT 
LB 

Kphl 


2.4011E+02 

7.0959E+01 

IN-LB 

Kpsi 

2.7692E-02 

-7 . 0959E+01 

2.4011E+02 

IN-LB 

Bz 



LB -SEC 

Bphi 


6 . 7528E-03 

-1.6341E-03 

IN-LB- SEC 

Bps! 


1.6341E-03 

6. 7528E-03 

IN-LB- SEC 
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( CASE 6 ) Face seal, no tilt, with pres, grad., lamda-10, sigma-200 
SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS - 1.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, 1.9099E+06 (RPM) 

INSIDE, OUTSIDE PRESSURE - 2.0000E+00, 1.0000E+00 (PS I) 

VISCOSITY- 8.3333E-11 (PSI-SEC), AMBIENT PRESSURE - 1.0000E+00 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 2 0 

CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES - -3.0782E-16, 1.1794E-15 (IN-LB) 

MINIMUM FILM THICKNESS - 1.0000E-03 (IN) 

FLOW - 2.2763E+01 (IN**3/SEC) MEASURED AT l.OOOOE+OO (PSI) 

TORQUE- 2.2645E-03 (IN-LB), FILM POWER LOSS - 6.8621E-03 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - 2.0000E+02 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

z (IN) 

2 . 38/2E+03 

phi (RAD) 

psi (RAD) 

FORCE UNIT 

Kz 

-4. 1732E-05 

-4. 1685E-05 

LB 

Kphi 

-2 . 5267E-08 

7 . 2912E+02 

-7.4833E+00 

IN-LB 

Kpsi 

5 . 8389E-08 

7.4833E+00 

7.2912E+02 

IN-LB 

Bz 

4. 1358E-03 

5 . 3850E-10 

5.3818E-10 
-5 . 8922E-05 

LB-SEC 

Bphi 

4. 2377E-13 

1 . 2426E-03 

IN-LB-SEC 

Bp si 

-5 . 5928E-13 

5 . 8922E-05 

1 . 2426E-03 

IN-LB-SEC 
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COMPARISON OF CASE 6 WITH PN471 


CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 

FLOW - 2 . 2763E+01 (IN**3/SEC) MEASURED AT 1.0000E+00 (PSI) 

TORQUE- 2 . 2644E-03 (IN-LB), FILM POWER LOSS - 6.8620E-03 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kz 

Kphi 

Kpsi 

Bz 

Bphi 

Bps! 


i$0E+< 


2.3870E+03 
4. 1383E-03 


phi (RAD) psi (RAD) 

7 . 2907E+02 -7.4987E+00 
7.4987E+00 7.2907E+02 

I.2438E-03 -5. 8939E-05 
5.8939E-05 1.2438E-03 


FORCE UNIT 
LB 

IN-LB 
IN- LB 
LB -SEC 
IN-LB-SEC 
IN-LB -SEC 


(HP) 
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( CASE 7 ) Misaligned shaft seal with pres. grad, to comp, w/ GASBEAR 
SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE - 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED - 1.9099E+05, 0.0000E+00 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 2.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY- 8 . 3333E-11 (PSI-SEC), AMBIENT PRESSURE - 1.0000E+00 (PSI) 

ITERATIONS AND ERROR CODE IN LAST PRESSURE CALCULATION - 3 0 

CALCULATED FORCES IN X.Y DIRECTIONS - 3.1535E+00. 1.2350E+00 (LB) 

CALCULATED MOMENTS ABOUt X,Y AXES - 3.7793E-01, -1.9883E-03 (IN-LB) 

MINIMUM FILM THICKNESS - 5.0000E-04 (IN) 

FLOW - 4 . 7315E+00 (IN**3/SEC) MEASURED AT 1.0000E+00 (PSI) 

TORQUE- 2 . 3248E-02 (IN-LB), FILM POWER LOSS - 7.0448E-02 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER - 0.0000E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 
x 


DISP. 

Kx 

Ky 

Kphi 

Kpsi 

Bx 

By 

Bphi 

Bps! 


9 . 22&7E+03 
-2 . 6167E+03 
-9 . 8594E+01 
7 . 1674E+02 


Y (IN) phi (RAD) psi (RAD) FORCE UNIT 

4.8183E+03 1.1084E+03 5.7870E+02 LB 

9 . 5964E+03 4.7033E+02 8.7124E+02 LB 

6 . 5068E+02 1.5189E+03 7.7760E+02 IN-LB 

-9 . 7210E+01 -1.0491E+03 1.2765E+03 IN-LB 

-2.0553E-01 -4 . 5816E-01 -1.0900E-01 6.2037E-02 LB-SEC 

2 . 5253E-01 -2 . 9852E-01 -6.7612E-02 -7.7454E-02 LB-SEC 

‘ 1.4706E-02 3.7737E-02 -8.3553E-02 IN-LB-SEC 

-5 . 8155E-02 5.7003E-02 6.7719E-02 IN-LB-SEC 


2 . 3869E-02 
-7.4893E-02 
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COMPARISON OF CASE 7 WITH GASBEAR USING SAME GRID AND ROMBERG EXTRAPOLATION 


CALCULATED FORCES IN X,Y DIRECTIONS - 3.1529E+00, 1.2362E+00 (LB) 

CALCULATED MOMENTS ABOUt X,Y AXES - 3.7751E-01, -7.3904E-04 (IN-LB) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Kx 



By 

Bphi 

Bps! 


FORCE UNIT 

LB 

LB 

IN-LB 

IN-LB 

LB-SEC 

LB-SEC 

IN-LB-SEC 

IN-LB-SEC 


3-57 



3.5 References 


1. Vohr, J.H. and Pan, C.H.T., 1 0n the Spiral Grooved Self Acting Gas Bearing', MTI-63TR52, 
Mechanical Technology Incorporated, Latham, NY, (1962) 

2. Vohr, J.H. and Pan, C.H.T., 'Design Data: Gas Lubricated Spin-Aids Bearings for Gyroscopes', 
MTI-68TR29, Mechanical Technology Incorporated, Latham, NY, (1968) 

3. Smalley, A.J., The Narrow Groove Theory of Spiral Grooved Gas Bearings: Development and 
Application of a Generalized Formulation for Numerical Solution', ASME J. Lub. Tech., V 94, 1, 
(1972), pp. 86-92 

4. Castelll, V. and Pirvics, J., 'Review of Methods in Gas Bearing Film Analysis'. Trans. ASME, 
(1968), pp. 777-792 

5. Press, W.H., Flannery, B.P., Teukolsky. S.A. and Vetterllng, W.T., 'Numerical Recipes', 
Cambridge University Press, (1986) 

6- Artiles, A., A., Walowit, J.A. and Shapiro, W., * 'Analysis of Hybrid Fluid Film Journal Bearings with 

Turbulence and Inertia Effects', Proc. ASME Symposium in Advances in Computer-Aided 
Design, (1982), pp. 25-52 

7 > Castelli, V., ‘Design of Gas Bearings - Volume 1, Part 4: Numerical Methods', Gas Bearing 

Course Notes, Mechanical Technology Incorporated, Latham, NY, (1971) 

8- Shapiro, W., ‘Computer Code SPIRALP for Gas-Lubricated Spiral-Groove Bearings and Seals', 

MTI 88TM2, Prepared for NASA under Contract No. NAS3-24645, (1988) 

9. Sato, Y., Ono, K. and Iwama, A., The Optimum Groove Geometry for Spiral Groove Viscous 
Pumps', ASME J. Tribology, V 112, 2, (1990), pp. 409-112 

1 0. Walowit, J. , * Users Manual For Computer Code SPIRALG Gas Lubricated Spiral Grooved Cylindrical 
And Face Seals', Mechanical Technology Inc. report 91TM11, September 1991 


3-58 



3.6 Nomenclature 


A dimensionless flow control area, (area/R 0 2 ) 

Aj coefficients of second order linear operator defined by Eq. (44), 1 = 1 ,6 
{a} column vector of eccentricity coefficients defined by Eq. (41) 
a k j kth component of {a} evaluated at grid point (i.j) 

B xy damping coefficient relating force in x direction to velocity in y direction, B^, B^, and B zz are 
similarly defined 

B^ damping coefficient relating moment about x axis to angular velocity about y axis, B^, B^ and 
B^^ are similarly defined 

B^ damping coefficient relating force in x direction to angular velocity about x axis, B^, B^, B^, 
B^ and B^ are similarly defined 

B *X damping coefficient relating moment about x axis to velocity in x direction, B +x , B^ y , B^ y , B^ and 
B^ z are similarly defined 

B^ dimensionless damping coefficient B^/Bq, same definitions apply to B^, B^, Byy, B zz 
dimensionless damping coefficient B^ + /(BgRg 2 ), same definitions apply to B^, B t(> , B^ t 
B^ dimensionless damping coefficient same definitions apply to B^, B^, B^, B^, 

B *X dimensionless damping coefficient B^ X /(B 0 R 0 ), same definitions apply to B +x , B^ y , B +y , B^, B^ z 

B 0 characteristic damping constant, I&RqVc 3 

[B] dimensionless damping coefficients in matrix form 

{b} column vector of coefficients of c' , arising from linearization of Eq. (39) 

C clearance (cylindrical seal) or reference film thickness (face seal) 

[d] coefficient matrix used in Newton-Raphson linearization procedure, see Eq. (34) 

[C 1 ] coefficient matrix obtained from steady state solution 
[d ,k ] derivative of [d] with respect to * k 

[d] diagonal coefficient matrix whose components are given by Eq. (51) 

[C* 1 ] complex coefficient matrix used in complex stiffness solution, [d] + &r[d] 

[d] coefficient matrix used in Newton-Raphson linearization procedure, see Eq. (34) 

[d] coefficient matrix obtained from steady state solution 
(d ,k ) derivative of [d] with respect to « k 

(E*) coefficient matrix used in Newton-Raphson linearization procedure, see Eq. (34) 

[E*] coefficient matrix obtained from steady state solution 
[E* ,k ] derivative of [E*] with respect to « k 
e x ,e y eccentricity in x,y direction (cylindrical seal) 
e z axial displacement, (face seal) 

Fj j residual outflow function from flow balance at grid point (i,j) 

G second order non-linear operator defined by Eq. (39) 

H steady state, unperturbed, value of H r 
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H r dimensionless film thickness, h r /C 

h g film thickness over grooves, see Fig. 2 

h r film thickness over ridges, see Fig. 2 

i,j,k subscripts used generically as indices 

i.l unit vectors in 0 ,s directions 
3 unit Imaginary number, /n 

stiffness coefficient relating force in x direction to displacement In y direction, K^, K^, Kyy and 
are similarly defined 

At stiffness coefficient relating moment about x axis to rotation about y axis, ■W A> and are 
similarly defined 

stiffness coefficient relating force in x direction to rotation about x axis, K^, K^, K^, and 
are similarly defined 

l^ x stiffness coefficient relating moment about x axis to displacement in x direction, K^ x , K^ y , K^y, 
K^ z and K^ z are similarly defined 

dimensionless stiffness coefficient K^/Kq, same definitions apply to K^, K^, Kyy, K^ 
dimensionless stiffness coefficient ^/(KqRq 2 ), same definitions apply to F^, K^, F^ 

K* dimensionless stiffness coefficient K^/OCqRq), same definitions apply to K^, K^, K^, K^, 

F^ x dimensionless stiffness coefficient ^/(KqRq), same definitions apply to F^ x , F^ y , F^ y , K^ r K^ z 

Kq characteristic stiffness constant, p 0 R 0 2 /C 

[K| dimensionless stiffnesses in matrix form 

k, spiral groove coefficient defined by Eq. (25), i-1,2_,8 

L seal length, see Fig. 1 

L dimensionless length, L/(2 Rq) 

l g groove width, rA6 g 

l r ridge width, rA0 r 

Sf second order linear operator defined by Eq. (44) 

M number of grid points In s direction 
M x ,M y applied moment about x,y axis 
tt x ,M y dimensionless moment, (M x .M y )/(R 0 3 Po) 

N number of grid points in 6 direction 
N g number of spiral grooves 

n unit vector normal to S 

np unit vector normal to groove 

P dimensionless pressure, (p-Pq)/p 0 

P steady state unperturbed value of P 

Pjj dimensionless pressure, P, at grid point (i,j) 

Pj dimensionless pressure, P, at point I shown In Fig. 3, 1=1 .....9 

P,,P f dimensionless boundary pressures (p l -p 0 )/p 0 , (p r -p^/p 0 
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{P* } column vector of dimensionless pressure disturbances due to perturbation in {* } 

{Pj } column vector of disturbance pressures associated with c k 
{P*} complex amplitude of {P' }, {P' } = {P*}e 3ot 
{Pj k } column vector of complex stiffness pressures associated with « k 
{P„ } real part of {P*} 

{P a } Imaginary part of {P*> 

p global pressure in absolute units 

P' local pressure In absolute units 

p 0 reference pressure in absolute units, normally taken to be the minimum of the boundary pressures 
p,.p r boundary pressures in absolute units at s,,s r 
6 dimensionless flow vector, I^RqCi/OjqC 3 ) 

Q S ,C^ components of dimensionless flow vector in sfi directions 
Qjj.Qjj dimensionless flow components shown in Fig. 3 

Q in dimensionless flow rate, I^^/^qC 3 ) 

q global flow vector, mass flow rate per unit transverse length divided by density at pressure p 0 

q s ,% components of q in s^ directions 

q A global mass flow rate per unit area displaced by rate of decrease of film, divided by density at p 0 
q in volumetric flow rate measured at pressure p 0 

O' local flow vector, mass flow rate per unit transverse length divided by density at pressure p 0 
q's.cjj components of q» in s,0 directions 

q A local mass flow rate per unit area displaced by rate of decrease of film, divided by density at p 0 
R dimensionless coordinate, r/R 0 , taken as 1 for cylindrical seal 

Rq reference radius, taken as outside radius for face seal and shaft radius for cylindrical seal 

Rj inside radius of face seal 

{R f } column vector used in Newton-Raphson linearization procedure, see Eq. (34) 

{R J > column vector obtained from steady state solution 

{Rl ,k } derivative of {R^} with respect to « k 

{Rj k } column vectors whose components are given by Eq. (52) 

{R j,k } complex column vectors used complex stiffness solution {R ,,k } -3*{R i,k } 
r radial coordinate, taken as R 0 for cylindrical seal 
S dimensionless coordinate, s/R Q 

S |( S r dimensionless boundary coordinates s,/Rq, s r /R 0 

S dimensionless length of line surrounding flow control area (length /R q) 

s transverse coordinate, s * r for a face seal and s » z for a cylindrical seal 

s g transverse coordinate at start of groove 

s, left boundary coordinate, s, = R, for face seal and s, * -L/2 for cylindrical seal 
s f right boundary coordinate, s r = R 0 for face seal and s r * L/2 for cylindrical seal 
T torque 
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T dimensionless torque, T/(p Q R 0 2 C) 

T c T obtained with course grid in Romberg extrapolation example 
T f T obtained with fine grid in Romberg extrapolation example 
T f T obtained by Romberg extrapolation 

t time 

1 dimensionless time, co t/(2A. ) 

tp unit vector tangential to groove 

u n velocity of grooved surface 

u 2 velocity of smooth surface 

W x ,W y applied loads in x,y direction (cylindrical seal) 

W x ,W y dimensionless loads (W x ,W y )/(p 0 R 0 2 ) 

W z applied load in z direction (face seal) 

W z dimensionless load, W z /(p 0 R 0 2 ) 

{W} column vector containing dimensionless loads and moments {W ,W ,M ,M } for cylindrical seal, 
{W z ,M x ,M y > for face seal 

{W} 0|d {W} at previous iteration in eccentricity homing process 
x,y coordinate variables, see Fig. 1 

Z dimensionless axial coordinate for cylindrical seal, z/Rg 

z axial coordinate, measured from axial center for cyl. seal or from reference film, C, for face seal 

a groove to pitch ratio, l g /(l g +l r ) 

a opt value of a for maximum stagnation pressure gradient in concentric cylindrical seal 

fi spiral groove angle, angle between grooves and surface velocity 

l opt value of fi for maximum stagnation pressure gradient in concentric cylindrical seal 

t numerical damping factor used In eccentricity homing process 
r fflm thickness ratio, h g /h r 

A A dimensionless flow control area about single grid point shaded area in Fig. 3 

A A, portion of A A in quadrant containing point i, see Fig. 3 

Ap global pressure difference, see Fig. 2 

A p» g pressure difference across groove, see Fig. 2 

A p t pressure difference across ridge, see Fig. 2 

A$|j line or arc length associated with Q|j 

ASjj line or arc length associated with Qjj 

A s transverse length of groove-ridge pair 

A8 circumferential extent of groove-ridge pair, see Fig. 2; (also used generally for change in 6) 

A6 g circumferential extent of groove, see Fig. 2 
A0 r circumferential extent of ridge, see Fig. 2 
i dimensionless groove depth, (h g -h f )/C 

£ opt value of J for maximum stagnation pressure gradient in concentric cylindrical seal 
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« x( f y eccentricity ratio, (e x ,e y )/C (cylindrical seal) 

< z axial displacement ratio e z /C, (face seal) 

« k kth component of eccentricity matrix 

(e] row matrix of eccentricity components, [« z *,t] (face seal) and [« x ,« y #,'H (cylindrical seal) 

{« } column vector, transpose of [c ] 

{< } old {< } at previous Iteration in homing process 
<' eccentricity disturbance function (scaler) 

ij small increment used in perturbations and in evaluating derivatives 
6 angular coordinate, see Fig. 1 

6 angular coordinate at start of groove 

A compressibility number, &<•> R Q 2 /{p Q C?) 

Aj groove compressibility number, A 2o>a (1 -o)sinn 

n viscosity 

a squeeze number, 2A a 

r global shear stress 

f dimensionless shear stress, tRq/(p 0 C) 

t' local shear stress 

* rotation about x axis 

$ reduced rotation, ♦ Rq/C 

7 rotation about y axis 

♦ reduced rotation 7 R 0 /C 

a angular velocity of disturbance 

a dimensionless disturbance frequency, o /<•> 

- <•> total angular velocity, to t + w 2 
a dimensionless angular velocity ratio, (u 2 -u ,)/« 

w 1 angular velocity of grooved surface 

2 angular velocity of smooth surface 

V gradient operator, (1/r)0/30)i + £/9s)], on dimensional quantities and (l/R)@/96)i + 
(3/9S)], on dimensionless quantities 
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4.0 Industrial Code ICYL-Incompressible, Cylindrical Seals 

Incompressible cylindrical seals are used to reduce leakage from higher pressures. 
The pressures generated in plain cylindrical seals with incompressible fluids 
typically result in forces which are normal to the displacement and therefore tend 
to destabilize the rotating shaft. Surface roughness, geometry alterations, and 
external pressurization are ways in which the direct stiffness and damping 
coefficients can be improved and the cross-coupled stiffness decreased in order 
to improve stability. 

The computer code ICYL was developed to evaluate the performance of cylindrical 
seals operating with incompressible fluids. The pressure and velocity distributions 
within the seal clearance are first evaluated from the governing equations. From 
these, design quantities such as seal leakage flows, power loss and resulting 
forces and moments are calculated. Minimum film thicknesses and maximum 
pressures as well as critical rotor-dynamics coefficients such as stiffness, damping 
and critical mass are evaluated. 

Program capabilities: 

1. 2-D incompressible isoviscous flow in cylindrical geometry. 

2. Rotation of both rotor and housing. 

3. Roughness of both rotor and housing. 

4. Arbitrary film thickness distribution, including features such as steps, 
pockets, tapers and preloaded arcs 

5. Rotor position described by four degrees of freedom (translational and 
rotational) 

6. Up to 32 dynamic coefficients as well as the critical mass may be calculated 
for use in rotor-dynamic design, including system response and stability 
calculations. 

7. External forces and moments may be prescribed independently to find rotor 
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position. 

8. Pocket pressures or orifice size are prescribed. 

9. Laminar or turbulent flow. 

10. Cavitation. 

11. Inertia pressure drop at inlets to fluid film (from ends of seal and from 
pressurized pockets). 

Assumptions 

1 . The film thickness is assumed to be small compared with seal lengths and 
diameters but large compared with surface roughness. 

2. Pockets supplied from an external pressure source through an orifice 
restriction are assumed to be sufficiently deep that the pressure is constant 
within them. 

3. Wall roughness is assumed to be isotropic and represented by an 
"equivalent sand roughness" height. 

4. Fluid inertia effects in the film are negligible. 
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4.1 THEORETICAL DESCRIPTION AND NUMERICAL METHODS 


Figures 4-1 through 4-3 illustrate the geometry of a cylindrical seal as well as the 
coordinate system used here to describe it. Figure 4-1 shows the seal housing of 
length L separated from the rotor by the film thickness H. The coordinate system 
is placed at the mid-length of the seal with the circumferential coordinate 0 
measured from the x-axis. Figure 4-2 illustrates the displaced, misaligned rotor, 
while figure 4-3 shows an axial cross-section of the film thickness. 


The film thickness and time rate-of-change thereof are written: 


H - H 0 - (e x +ZB) C 0 #e - (« y -Z4)sin0 


dH 

dt 


loose - 



Sine 


( 4 - 1 ) 


where H 0 , an arbitrary function of film coordinates (B,Z), represents the film 
thickness distribution for a rotor that is aligned and centered with the housing. e x 
and e y represent the components of rotor eccentricity at the seal mid-length, while 
A and B represent the angles of rotor rotation about the x and y axes, respectively. 
The former are referred to as the radial or lateral displacements and the latter as 
the angular displacements. 


Governing equations 

The equations governing the flow of incompressible fluids in thin films are obtained 
[10,11,15] by integrating the Navier-Stokes momentum and continuity equations 


4-3 


CYLINDRICAL SEAL GEOMETRY SCHEMATIC 

(CONCENTRIC ALIGNED POSITION) 



Figure 4-1 Cylindrical Seal Geometry 
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SCHEMATIC ILLUSTRAT 


Figure 4 





Rotor with Lateral and Angular Displacements 
i-5 


AXIAL CROSS SECTION OF SEAL 
WITH ECCENTRIC ROTOR 

(FILM THICKNESS EXAGGERATED) 



Figure 4-3 Axial Cross-Section of Seal with Eccentric Rotor 
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across the film 1 : 


( fjRe^Re^ [J h* qP + {Re^U^ Re b f b U^ 

2 " ~ \iR ae 2 

(4-2) 

(fjRej+ftReJ ^ _tf_dP 
2 " ~ |x dZ 

)f£m * * f - ° <«) 

where fj and f b are the friction factors relative to the housing and journal surfaces, 
respectively, and are functions of the Reynolds numbers relative to these surfaces 
as well as of their roughness. They are given by: 

Re,- 2-^J(U-U,) z +V* (4-4) 

v> 


where i=j,b, and: 




~ , Re,i 1000 (lamina/) 
-3S l +2£*) + / r / *(3{ 2 -2{ 3 ), 1000</?e / <3000 

nVy 

t; , Re, i 3000 (turbulent) 


(4-5) 


( ■ 


Re ,- 1000 
2000 


0.001375 


1 a. 

'lO 4 *, 10* ' 

jf 

3 

1 ▼ 

H + 4/teJ 



(4-6) 


The friction factor for turbulent flow through pipes, f*, in equation (4-6) uses the 
curve-fit obtained by Nelson [12] to Moody’s data. The transition from laminar to 


1 the word film or the term film thickness will be used to mean the gap of lubricant separating the 
rotor and housing. 
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turbulent flow is obtained using a cubic polynomial which matches values and 
slopes at both ends, as reflected by equation (4-5). Figure 4-4 is a plot of the 
friction factor versus Reynolds number and surface roughness, while Figure 4-5 is 
an enlargement showing the detail of the transition region. 


Under laminar flow with the friction factors equal to 12/Re, the velocities can be 
solved explicitly in terms of the pressure gradients: 

u- -12 + A/ + A b 

de 1 b 

( 4 - 7 ) 

v= -12 h 2 & 

dz 


Lubrication Background: 

In the classical theory of lubrication, when the housing is stationary and the rotor 
wall velocity is Uj=oR, the fluid velocity components are expressed explicitly in 
terms of the pressure gradients: 

Um + Vm («) 
12)i A 86 2 ’ 12|i dz 


where G x and G 2 are turbulence coefficients^] which become unity in the laminar 
regime. Substituting these velocity components into the continuity equation, 
results in the classical Reynold’s equation: 




Boundary conditions 

Boundary conditions on the film pressure distribution consist on prescribing either 
the pressure at the boundaries of the film, the flow normal to these boundaries, or 
a relation between these two quantities. 
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Selected friction factor (Moody) 



Reynolds number 


Figure 4-4 Friction factor versus Reynolds number 
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Friction (odor, f 


Hg.1 Transition friction factor 



Re = pVh/|t 


Figure 4-5 Detail of friction factor in transition region 
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At the circumferential ends of the seal surface model, either the pressures are 
prescribed: P( Z, 0J = 0 and P( Z, 0,)= 0, 

or periodic boundary conditions exists: 

P(Z,0 8 )= P(Z,0.) and U(Z,© 3 )= U(Z,0 # ). 

Periodic boundary conditions are used, for example, for a 360* seal, where 
0,=0,+2ft. 

At the left end of the seal surface model, the pressure/flow relationship is 
prescribed: 

P(-L/2,0). P,-K,! ! pV„ 2 . 

At the right end either the same relationship is used: 

P(L/2,6)- P r -K.>!pV„ 2 , 
or the axial velocity is set to zero: 

V(O,0)=O 

when a symmetry boundary is present at the seal mid-length. 

Finally, at pocket boundaries: 

P(Z,8)= P p - K, VpV„ 2 . 


In all of the above relationships, 

P-A P/»>0 
V„ • 

l 0, 9-AiO (4-10) 

P ■ Ui z + Vi, 

is the flow velocity at the entrance to the film, normal to the pressurized boundary. 
No pressure drop exists in the case of reverse flow (i.e., flow into the pressurized 
boundary). 
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External pressurization 

The pressure drop across the orifice supplying the pocket is given by: 


>f(^r ,4 - ,i) 

where A 0 is the orifice area, C d is the discharge coefficient and the flow Q r is 
obtained by satisfying continuity over the pocket volume: 

Q r -fH9AdS + f2gdA (4-12) 

i, J A, dt 

where Ap is the pocket area, S p is its perimeter. Note that the contribution of V« ft 
to this last equation may be positive or negative. 


Dimensionless variables 

Using the following transformation to dimensionless variables, 


b- B (C 3 /12 mR 4 ) 
>- F/(P 0 R 2 ) 
h= H/C 


A„- 6pU B R/(C 2 P 0 ) 
A,- BmU.RAC^PJ 
€ = e/C 


k= K (C/P 0 R 2 ) a - A (R/C) 

m- M/ (P 0 R 3 ) 0 = B (R/C) 


P= P/P 0 

q I -Q r (12p/P 0 C 3 ) 
u= U (12 hR/C 2 P 0 ) 
v= V $2 hR/C?PJ 
z = Z/R 

r - t (C 2 P 0 /12 fiB 2 ) 


Re'- phVp/p 2 
Re 0 '= pC^/fR,* 2 ) 

A,= pC*P 0 /(288 A 0 ?C d 2 p 2 ) 

= (Re 0 ’/288) (C 3 R/A 0 2 C d 2 ) 
A.- X. (Re„'C/288 R), 
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equations (4-1), (4-2), (4-3) and (4-4) become: 

h » h 0 - (e x +z(})cos0 - (e r -z«)sin6 

^z^-lcose - jslne 

dx) dx) 


dh 

dx 


dt. 


dx 


( 4 - 13 ) 


- 12 /,* J| * {Re/fi. Re t f t AJ 


VlROi+f'.ReJ i£ 

2 az 


( 4 - 14 ) 


i<“"» * S w f ■ 0 <*■'*> 

Da * h 

Re,- (u- 2A,) ! ♦ i* / - /./> (<-«) 

Equations (4-5) and (4-6) remained unaltered, as they were already dimensionless. 


The dimensionless form of the boundary conditions now become: 
At the circumferential ends, either: 
p(z, ej= 0 and p(z, 0 # )= 0 
or: 

P( z, 0.) = P( z, 6.) and u( z, e,)= u( z, 6 # ). 
when periodic boundary conditions are present. 

At the left end: 

P(-L/D.8)= P, -A.v„ 2 , 
and at the right end either 

p( L/D,6)= p, -A.»„ 2 
or: 

v(O,0)=O. 
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At pocket boundaries: p(z,0) = p p - A # v n 2 

where: 

fp-A 9 A>0 

V * "1 0 , ( 4 - 17 ) 

V * u$ M +vi 9 

Equations (4-11) and (4-12) governing the external pressurization become: 

P,~ P p m sgrKq,)A f q? ( 4 - 18 ) 

q, - jh v h ds + dQdz ( 4 - 19 ) 

• A ^ 


Solution of film pressures 

Discretization of the seal surface is done by using a rectangular grid, with M lines 
in the axial direction and N lines in the circumferential direction. The grid lines are 
separated by variable increments. The pressure distribution is represented by 
discrete values at the grid points located at the intersections of the grid lines. 
There must be grid lines coincident with the boundaries of the seal surface 
(Z=± L/2, 0 =0,, 0 =0,) and with the pocket boundaries. Using the cell method [3], 
a control area or cell is centered at each grid point and extending half way to the 
neighboring grid lines, as shown by the shaded area in Figure 4-6. The grid points 
are noted by the solid circles and have grid coordinates i.j. The film thickness is 
evaluated at the corners of the cell (denoted by the shaded circles marked h,, h 2 , 
h 3 , and hj located at the geometric centers of the rectangles formed by the grid 
lines. This staggered configuration allows a discontinuous film thickness to be 
treated, as occurs, for example in a seal with a Rayleigh-step. Circumferential and 
axial components of velocity are also associated with each of the four cell corners. 

Using the divergence theorem, the continuity equation may be integrated over the 
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Figure 4-6 Flow control area about grid point i,j. 
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cell to give: 


-f vh dS - f ^ dA ( 4 . 20 ) 

a, a, 

where \ and S c are the cell area and perimeters, respectively. The left hand side 
of the above equation is the sum of the flows out of the cell while the right hand 
side is the rate of change of the cell volume. The finite-difference form of this 
equation is: 

F V - -y ! ( 0 i A i ~ u * h 4 ) ♦ -Mb) ♦ 

♦ -MW + -v*>h) - <«<> 

' 4 J? (**i***nt*«l***n) • <> 

where F (J is the error in satisfying continuity of flow in the cell centered at i.j. 
Although the time rate of change of film thickness has been evaluated at the center 
of the cell, it could have alternatively been evaluated at each of the four cell 
comers. 


When the grid point falls on a pressurized boundary, such as a pocket or seal end, 
the film pressure error is: 


F u m Pb~ Pu -A # max(0 l v ffl ) 2 - 0 



(4-22) 


where p b is the dimensionless boundary pressure 2 , v n is the mean velocity of the 
flow that crosses the portion of the boundary perimeter that intersects the cell, and 
£|j represents the sum of the appropriate terms in equation (4-21) contributing to 
the cell flow. Figure 4-7 shows an example of the cell i.j located at the right bottom 


2 
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corner of a pocket. In this case, the mean velocity would be evaluated as: 


v « - 


Az, 


(u,h,) * -u,h,) * v,h , -KjA,) - ^-(v 3 h,) - 


A0, 


dhy ( A Z/ + A Z M ) A 8 y + A Z M A By^ 

0T 4 


A 

(A8 y . 1+ A z,)h u 

▼ 

2 J 


(4-23) 


Equations (4-21) and (4-22) represent the finite-difference form of the continuity 
equation that must be solved for the pressures. The eight components of velocity 
used in these equations are functions of the nine pressures at or neighboring grid 
point i,j, and are evaluated as described in section 2.3. Following the procedure 
described in reference 1, these highly nonlinear equations can be solved using the 
Newton-Raphson iteration method [14]. The procedure is started with an initially 
guessed or previously calculated pressure distribution, p (J . The error function F § 
is then linearized about this guess in order to obtain a better approximation to the 
pressures P,/’"": 

F, + E -^(Pir-Pj-o (4-24) 

I-I-1.H 

where a forward difference or a central difference may optionally be used to 
numerically evaluate the partial derivatives. Pressures without the superscript new 
relate to the previous or "old" approximation. If we introduce the column vector 
{pj"*"} as the M new pressures at the Jth column of grid points, Equation (4-24) 
may be written: 

Iduri ♦ IE'Hr-T) * [D'HP |T) - {«'} - (4 " 25) 


where |C J ], [E J ] and |D*] are tri-diagonal matrices whose interior elements are: 




^1*4 






. Dj* 


^Pl**J* 


k- -1,0,1 ; I -2,...,M-1 . 


The interior elements of the column vector {R 1 } are: 
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1 

R|* ■ £ (^U-kP|»k,J + E|!l*kP|«k,J-1 + D|!l»kPl.k.J»l) “ F|j • 

k --1 


The set of linear equations (4-26) that result for the next guess of pressure 
distribution is in a form suitable for solution by the column method which is 
described in detail in References 3 and 4. This method makes use of the banded 
nature of the equations in order to minimize computer time. 


Solution of flow velocity 


The momentum equations (4-14) are used in order to evaluate the velocity 
components from the pressure gradients. These equations may be rewritten in the 
generic form: 



!l^ds5h u * 12 />*|f - (Re/faRe^ A.) - 0. 

* ( 4 . 26 ) 

f j Re i +f b Re b v + 12 /,2i£ , o, 

2 dz 


where the Reynolds numbers used to evaluate the friction factors are based on the 
magnitude of the local fluid velocity relative to each surface: 


R»r ("- 2 A/ * V s , 

( 4 - 27 ) 

( U ~ 2A bf + ^ 2 . 


The dependence of the friction factors on velocity components orthogonal to each 
momentum direction couples the two momentum equations. Figure 4-8 is a 
schematic of the rectangular region between axial grid lines i and i+1 and 
circumferential grid lines j and j+1. In order to preserve continuity, it is essential 
that the same equation be used to evaluate the velocity components for adjacent 
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Figure 4-8 Schematic of rectangular region between grid lines. 


4-20 






cells. That is, the velocity u, out of the shaded cell centered at i,j must have the 
same value as the velocity u 4 into the cell centered at i,j+1. This value is 
designated as u' in the figure. Similarly, the velocity v, out of the cell i,j must be 
the same as v 2 into the cell at i + 1 ,j, and is designated as v*. This is achieved by 
using the average of the two corresponding orthogonal components. Thus, the 
component u* is determined by the u-momentum equation: 



Pij+\ ~Pu 


u~. 


v~ + 
2 


= 0 


while the component v’ is determined by the v-momentum equation: 


( 4 - 28 ) 




Az, 


, v 


= 0 


( 4 - 29 ) 


Similarly, u + and v + are determined by: 


Pmj^~Pi*aj v~ + v* 
• ~ » u » z 


A 6 



Pi*\j+\~Pij+\ 
A z i 


u'+u* 




= 0 


( 4 - 30 ) 


Equations (4-28), (4-29) and (4-30) are four coupled equations that determine the 
velocity components from the four pressures at the corners of the rectangle 
between grid lines and must be solved simultaneously. This is done using an inner 
Newton-Raphson iteration loop. By performing the differentiation of the error 
functions (Gu, Gv, ...) with respect to the four unknown velocities, analytically 
instead of numerically, significant computer time is saved. If the velocities have not 
been previously calculated initial guesses may be obtained from equations (4-7) 
assuming laminar flow. Once the iterations for the velocities have converged, their 
values are saved to provide a good starting guess for the next time they must be 
calculated. 


One simplification is possible by assuming that the friction factors are constant 
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within the rectangular region and the Reynolds numbers are based on the 
averaged flow velocity components, h (u'+u + ) and h (v +v + ). Although this does 
not uncouple the four equations, it requires less number of evaluations of the 
square root in equation (4-16). Since this assumption saves some computer time 
without introducing significant errors, it was chosen as the default program option 
(IFRIC = 3). However, occasionally when the grid is not very fine and the pressure 
gradients vary rapidly, the iterations will diverge and the more rigorous formulation, 
which uses distinct friction factors for each of the four momentum equations, 
should be used with the IFRIC =4 option. 

If the surfaces are smooth and the housing is stationary so that the continuity 
equation takes the form of equation (4-9), the simpler formulation described in 
detail in Reference 1 may be used by selecting the option IFRIC =0, resulting in 
significant reduction in computer time. 


Fluid film load, moment and torque 

The forces and moments on the rotor generated by the fluid film pressure 
distribution are obtained by integration of the pressure distribution over the 
cylindrical seal surface: 


\ F . 

F y 

L •• 

0090 

sine 


M x 

. // p 

-L 6, 

-ZsinO 

Rdd dZ 

My 


ZoosO 



The dimensionless form this equation is written: 


for] 

L 

0090 


f y 

D •• 

sine 


* 

m x 

- / / p 

L 6. 

-zsine 

dQ dz 

(fly 

D 

zoosO 



( 4 - 31 ) 


( 4 - 32 ) 


The differential of torque transmitted from the housing to the rotor is given by the 
cross product of the position vector T and the shear traction vector acting on the 
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housing 1: 

For laminar regime, fjFtej = f b Re b = 12, and the equation simplifies to: 


f = Te z = / J> x fdA 

A t 

’ R jf S r x (*•*,) dA 

T ■ ffl*S - <*cz 


( 4 - 33 ) 


P n R* 


T = 

2C e i/' ae 


//<*£ 




( 4 - 34 ) 


The power loss due to the difference in velocities across the two surfaces is 
obtained by doting this torque with the relative velocity: 


P 0 R* 


2C. 


P « T(v> b - «y) 

//|*S - (VAt) ^ 


( 4 - 35 ) 


O A. 


72 h 


Stiffness and damping coefficients 

Defining W to be a generalized vector of forces and moments generated by the 
fluid film pressure and T to be a generalized vector of lateral and angular 
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displacements: 


f r 

r* 

«* 

e y 

• f-i' 


m x 


a 

dx 

a 

m y 


P. 


Ip) 


the matrices of stiffness and damping coefficients can be written: 



aw, 

dt, 


( 4 - 36 ) 


( 4 - 37 ) 


where the subscripts i and j range over x, y, a and p. These coefficients are 
evaluated by numerical differentiation of W, using a forward difference. For 
example: 


K V P) ~ ^( c x» c » P) 

* * ft 


( 4 - 38 ) 


Solution of rotor position and pocket pressures 

If the rotor position is specified, equation (4-36) is used to solve for the fluid film 
forces and moments in terms of the calculated pressure field. Similarly, if the 
pocket pressures are specified, equation (4-1 1) is used to solve for the orifice size 
in terms of the supply pressure and calculated pocket flow. 

On the other hand, if externally applied loads and moments on the rotor (f,^, 
and m^ are specified they must be balanced by the fluid film forces to 
maintain static equilibrium. Similarly, once the orifice size is specified, equation 
(4-11) must be satisfied by the pressure in each pocket. The global set of 
equations that must be satisfied by the rotor displacements and pocket pressures 
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are: 


m - 
m - - f » 

rnjA - ~ m xg 
m^i) * -rrtyg 

P, - P p1 - sng^Pn) Ki (<7„f . for pocket 1 , 

P,~ Ppt m for pocket 2, etc. 


( 4 - 39 ) 


The vector T can now be redefined to include the pocket pressures and a 
generalized vector of errors in forces, moments and pocket pressures W # can be 
defined: 


r ' 

€ v 



X 

r 


f y + f » 

z y 

a 


ff^x ^ ff^xg 

P 

w 9 - 

m y + m n 

Ppi 


P. - Ppi - sC'Kffrt) 

Pp2 


P,-Ppi- sgfi(<l^*a(<leT 


( 4 - 40 ) 


Solution of the global equations is performed by Newton-Raphson iterations, as 
follows: 

( r*"" - r ) « 0 I 4-41 ) 

where, as before, the superscript new indicates the newer values of the vector r. 


K 


dW, 


dr 
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4.2 SAMPLE PROBLEMS 


A number of sample problems have been prepared to demonstrate the behavior 
and various features of the computer program. They are intended primarily for 
illustration and do not necessarily represent recommended seal designs. Table 4-1 
summarizes the mesh size, approximate execution times (on a 20Mz 386 PC) as 
well as a list of what variables where specified and iterated for in the outer loop for 
the sample problems. 

Samples EX1, EX2 and EX3 were selected with a coarse (5x11) mesh covering a 
90* sector in order to demonstrate the use of pressurized pockets and iterations 
for rotor position within a reasonable execution time. A pocket with a supply 
pressure of 100 psi was centered on the seal sector modeled. 

Sample EX1 contained two cases. In the first case, the pocket pressure 
was specified as 50 psi, resulting in an orifice diameter of 0.0137 inches 
calculated by the program. Both components of the resulting fluid film force 
are equal and the moments are zero, as would be expected at the 
concentric position. In the second case, the rotor was moved with to 
eccentricity ratio of e x = 0.1 and given a misalignment ratio of/? =0.1 about 
the y-axis. With the value of orifice diameter already assigned from the first 
case, the pocket pressure and forces rise slightly, generating non-zero 
moments. 

In sample EX2, external forces and moments equal to the negative of those 
resulting in EX1 were specified, in order to have the program iterate for the 
rotor radial and angular positions. Five unknowns, the four displacement 
components (e x ,e y , a, and /9) as well as the pocket pressure, are iterated 
for simultaneously. Although it wasn’t needed, IREADP = 1 and READP = 
’ICYLEX1.888’ where specified in order to illustrate the use of reading the 
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Table 4-1 


Summary of sample cases 


Case Mesh ISYM Variable App.exec Features 



Size 


Found’ 

Specified 

time” 


EX1 

5x11 

0 

dorif .Pport 

Ppocn.E«»dorif 

4.6 min 

1-pocket 

EX2 

5x11 

0 

P.Ppoo. 

F«,F yi M #t My,dorif 

11.6 min 

1-pocket 

EX3 

5x11 

0 


M,,My f dorif 

6.2 min 

Tapered pocket 

F3 

9x61 

1 

- 

all 

29 min 

Raleigh-step 

F4 

7x61 

0 

- 

all 

7.8 min 

Axial taper 

F9 

5x73 

1 

K,B 

all (3 preloads) 

1.6 hrs 

3-lobe 

II 

5x61 

1 

dor if 

P^ 

7 min 

4-pocket 

12 

5x61 

1 


£,,dor1f 

1.8 hrs 

4-pocket 

13 

5x61 

1 

e., 

F*,F,,dorif 

1.9 hrs 

4-pocket 

14 

9x61 

0 

K.B.P^c 

€, t B f dor1f 

7.7 hrs 

4-pocket 

15 

11x61 

0 

dor if 

P— 

5.2 min 

8-pocket 

16 

11x61 

0 

P— 

dor if 

3.1 hrs 

8-pocket 

015 

5x31 

0 

K.B 

all 

1 hr 45 sec 

Roughness 


Indicate whether stiffness, damping coefficients were requested, 
on an IBM PS/2 model 70 (386 20-Mhz) computer. 


pressure distribution from a previous run. Since the iteration was begun at 
the concentric position where the orifice was sized, the pocket flow error 
was zero and increased when the rotor was moved in the first iteration, 
causing the run to abort. When the limit on diverging iterations (MAXDIT) 
was increased to 2, the iterations converged in only 3 iterations to within a 
small error of the values expected (e x = 0.1, p =0.1). 

In sample EX3, an axial taper of ±30% of the clearance from end to end 
was superimposed. This calculation might be desirable to see the effect of 
machining tolerances or imperfection on seal components. This was 
accomplished by decreasing the clearance by 0.0003 inches 
(DELTA(1,1) =-0.003) as well as using DELTA(2,1) =0.006. In this sample, 
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the program was asked to find the angular rotor position such that no 
external moments where required, while the rotor eccentricity was varied, 
using e x = 0.1, 0.3 and 0.5 for the 1st, 2nd and 3rd cases, respectively. 
The results show that as the eccentricity is increased in the x-direction, the 
rotor twists about the y-axis in order for the moments generated by the film 
to be zero, resulting in >9 =0.050 and 0.17 at e x = 0.1 and 0.3, respectively. 
The case of e x =0.5 resulted in a negative film thickness with the appropriate 
error message and recommended user action: 

o REDUCE THE SPECIFIED APPLIED FORCES/MOMENTS 
o REDUCE THE SPECIFIED ECCENTRICITY/MISALIGNMENT 

The resulting film thickness is shown in Figure 4-9 and the pressure 
distributions are shown in Figure 4-10. 

Sample F3 shows a 120* sector with a Raleigh step of linearly varying 
depth. The resulting film thickness and pressure distributions are shown in 
Figures 4-11 and 4-12, respectively. 

Sample F4 shows a 1 20* sector with an axial taper in the right half (4 < i < 7) 
of 0.001 inches. Since two less intervals were used in the axial direction 
than in the previous cases, and since half as many iterations were required 
for the pressure distributions, the execution time was reduced from about 
29 to 8 minutes. The film thickness and pressure distributions are shown 
in Figures 4-13 and 4-14, respectively. 

Sample F9 is that of a full 360* seal with three 60° lobes. The dynamic 
coefficients was requested as the preload was increased from 0.1 to 0.3 in 
the middle case, to 0.8 at the last case. The film thickness and pressure 
distributions are shown in Figures 4-15 and 4-16, respectively. 
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ICYL: ICYLEX3 FILM TH I CHNI 
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Figure 4-9 Film thickness distribution for sample 
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Figure 4-10 Pressure distribution for sample 
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Figure 4-11 Film thickness distribution for sample F3 
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Figure 4-12 Pressure distribution for sample F3 
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Figure 4-14 Pressure distribution for sample F4 



ICYL: ICYLF9 FILM THICHNESS DISTRIBUTION 
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Figure 4-15 Film thickness distribution for sample F9 
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Figure 4-16 Pressure distribution for sample F9 
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Figure 4-17 Pressure distribution for sample II 




Samples II through 14 are more realistic models since they represent the full 
circumference. In sample II, the program takes about 7 minutes to 
calculate the orifice diameter for given pocket pressures. The resulting 
pressure distribution is shown in Figures 4-17. Sample 12 is the same as 
II except that the eccentricity and orifice diameter were prescribed, 
requiring the program to solve for the four pocket pressures. The resulting 
film thickness and pressure distributions are shown in Figures 4-18 and 
4-19, respectively. In sample 13, the orifice diameter as well as the radial 
force were prescribed, requiring the program to solve for the radial position 
as well as the pocket pressures. 

Sample 14 shows the dramatic increase in execution time with the number 
of axial grid lines, M. Sample II is a model of only half of the seal 
(ISYM = 1) at the concentric position with the pocket pressures specified. 
This run executes in less than 7 minutes in spite of the 5x61 mesh. 
Samples 14 is a model of the full axial length (ISYM=0), with an 11x61 
mesh, in which non-zero c x , a and orifice size are specified and all 32 
dynamic coefficients are requested. This run took 7.7 hours to execute. 
Calculations for each of these coefficients require convergence of the outer 
iteration loop with four unknown pocket pressures. 

Sample 15 and 16 are models representing the full circumference and length 
with two rows of 4 pockets. The orifice size is calculated in the concentric 
aligned position in 15 while 16 calculates what happens when the rotor is 
displaced to the e x = 0.4 position and rotated about the x-axis by a =0.4. For 
15, the resulting pressure distribution is shown in Figure 4-20. For 16, the film 
thickness and pressure distributions are shown in Figures 4-21 and 4-22, 
respectively. 
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Figure 4-19 Pressure distribution for sample 12 
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Figure 4-20 Pressure distribution for sample 15 
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Figure 4-21 Film thickness distribution for sample 16 
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Figure 4-22 Pressure distribution for sample 16 



Critical Mass (Lb-sec ^in) 


Sample 015 is a case of a plain cylindrical seal with increasing housing wall 
roughness. The wall roughness (ROUGHB) was varied from 1x1 O' 6 to 
IxlO' 3 inches in a logarithmic scale using IPAR = 14 and NPAR=-4. This 
input was used to generate the top curve of critical mass versus roughness 
shown in Figure 4-23. The stabilizing effect of housing roughness is more 
pronounced at the higher pressures due to the increased effect of inlet 
inertia. The last lines in the input file show how one would run additional 
values of pressure within the same run. These were not run because they 
were placed after the line with IST0P=1 in order to reduce the size of the 
output file. 



Housing roughness e/C 

Figure 4-23 Critical mass versus housing roughness 
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4.3 VERIFICATION 


ICYL has been compared with the results of two other MTI computer codes as well 
as currently published data. The first comparison was against a generic bearing 
program with many similar capabilities (GBEAR) based on the turbulent lubrication 
theory of Ng and Pan. A second comparison against a laminar bearing program 
(GASBEAR) was used to verify the calculations of moments and angular 
coefficients. Finally, comparison were made against calculations published by San 
Andres in Reference [16]. 

Comparison against MTI other codes 

The first of the MTI computer codes is GBEAR which is fully described in 
Reference 1. This program is based on the turbulent lubrication theory of Ng and 
Pan [13], and does not include surface roughness, housing rotation or calculation 
of misalignment coefficients. It includes inertia pressure drop at exit from pockets 
but not from the seal ends. 

Calculations were made with a 90* seal sector at an eccentricity ratio of 0.5 and 
with a pocket at its center with a prescribed pressure ratio of 0.5. Table 4-2 shows 
a comparison of pocket flow, orifices size, force, and stiffness and damping 
coefficients. As expected, comparisons of GBEAR against ICYL with the same 
friction model (IFRIC=0) yielded nearly identical results. With the new friction 
model that includes surface roughness effects, ICYL calculates lower torque(-32%), 
lower pocket flow (-13%) and orifice size (-7%), and force components (-6%). Very 
good agreement in the stiffness coefficients (-4%), and slightly higher damping 
coefficients^ 13%) are obtained. 

Other comparisons against GBEAR in the laminar regime and without pockets 
yielded identical results. 
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Table 4-2 Comparison against GBEAR. 



GBEAR 

ICYL 

IFRIC=0 

ICYL 

IFRIC=3 

ICYL 

IFRIC=4 

Recess flow (in 3 /s) 

25.75 

25.21 

20.931 

22.316 


0.0833 

0.0820 

0.0752 

0.0776 


14.38 

14.32 

8.791 

9.771 


1 


27.617 

brim 

Fx (Lb) 

3,694 


BRIM 

HHB9 

H (Lb) 

-3 , 488 

-3, 122 



Kxx (10 s Lb/ in) 

2.352 

2.267 

2.329 

2.344 


-1.461 

-1.378 

-1.280 

-1.397 


-1.998 

-1.874 

-1.871 

-1.961 


1.573 

1.481 

1.406 

1.564 

Bxx (Lb/in) 

232.08 

234.79 

269.01 

274.46 

Bxy (Lb/in) 

-175.53 

-175.87 

-194.38 

-199.65 

Byx (Lb/in) 

-174.78 

-174.10 

-192.40 

-200.56 

-®yy (Lb/in) 

173.87 

173.79 

187.57 

196.53 


A second MTI computer code with the fluid compressibility turned off (GASBEAR) 
was used to verify the calculation of the 24 stiffness and damping coefficients 
which involve rotor misalignment. GASBEAR was written for use in conjunction 
with plane journal bearings and cylindrical seals and does not treat turbulence or 
pressurized pockets. The comparison, in the laminar regime and with the same 
finite difference mesh, yielded identical coefficients. 

Comparison against published data 

A detailed comparison was made of the 5-pad hydrostatic bearing discussed by 
San Andres in Reference [16]. This high speed hybrid journal bearing operates 
at relatively high levels of pressurization and relatively low viscosity lubricants, in 
which the effects of pressure-induced turbulence become important. Fluid inertia 
may also be important. Figure 4-24 is a plot of the pressure distribution at the 
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concentric position, while Figure 4*25 and Figure 4-26 plot it for 40% eccentricity 
ratio of the journal between pockets and over a pocket center, respectively. 
Reproductions of the corresponding pressure distributions published by San 
Andres are included in the figures for comparison. It is noticed that the size of the 
pressure drops at the pocket exits (i.e., entrance to the film) as well as the general 
pressure distribution are comparable for both analyses. 

At the concentric position, bearing flow requirements calculated by ICYL is 42 
versus about 44 l/min reported by San Andres. Figure 4-27 and Figure 4-28 are 
plots comparing the direct and cross coupled stiffness coefficients while, 
Figure 4-29 and Figure 4-30 compare the direct and cross coupled damping 
coefficients, respectively, versus eccentricity ratio. In general, ICYL predicts about 
35% higher direct stiffness, 10% lower cross coupled stiffness coefficients, and 
15% lower direct damping at the concentric position. With increasing eccentricity 
ratio, the coefficients are observed to behave similarly and some of the 
discrepancies decrease. The cross-coupled damping coefficients with ICYL are 
equal in magnitude, opposite in sign and zero at the concentric position, as is 
expected with an incompressible fluid. San Andres’ non-zero concentric value (60 
kN-s/m) is due to the fluid compressibility in the pocket. Figure 4-31 shows the 
critical mass versus eccentricity. The concentric value of 119 Kg shows better 
stability than predicted by San Andres, which predicts an unstable bearing with a 
30 Kg mass. 

The analysis of San Andres includes the effect of fluid inertia in the film as well as 
some special effects inside the pocket, such as fluid compressibility and a one- 
dimensional circumferential pressure rise downstream of the orifice. There is also 
a slight difference in friction law used: MTI’s analysis follows the formula derived 
by Nelson[12] for Moody diagram, in which the term containing the Reynolds 
number is raised to the 1/3 power while San Andres uses the same formula with 
the power changed to 1/2.65 for a more restricted range of Reynolds numbers. 
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Figure 4-24 Comparison to SanAndres' 5-pad bearing at concentric position 
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Figure 4-25 Comparison to SanAndres’ 5-pad bearing with 40% eccentricity 
between pockets 
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Figure 4-26 Comparison to SanAndres’ 5-pad bearing with 40% eccentricity over 
pocket 
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Comparison of direct stiffnesses 
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Figure 4-27 Comparison of direct stiffness coefficient 
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Figure 4-28 Comparison of cross stiffness coefficients 



4-52 



Comparison of direct damping 
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Figure 4-29 Comparison of direct damping coefficients 
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Comparison of cross damping 



Eccentricity ratio Ex 


Figure 4-30 Comparison of cross damping coefficients 
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The above comparisons should provide reasonable verification, as the only 
discrepancies between the results can be explained by the different friction models 
and features between the codes. 
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4.5 NOMENCLATURE 


A, B 

b fj = By (C 3 /1^R 4 ) 
C 

e b- e j 


misalignment of rotor about the x and y axes, 
respectively, [radians] 

dimensionless damping coefficient matrix, where 
i,j= x, y, a, (3. 

nominal clearance. [L] 

roughness of the housing and journal surfaces. 

[L] 


f = F/(P 0 R 2 ) 
H 


components of rotor eccentricity at Z=0. [L] 

rotor eccentricity at Z=0. [L] 

components of fluid film force about x and y 
axes. [F] 

dimensionless fluid film force, 
local film thickness. [L] 


H 0 

h = H/C 


local film thickness for the concentric aligned 
rotor (i.e., e x = e y = A = B = 0). [L] 
dimensionless local film thickness. 


K. 


coefficient of pressure drop at inlet to film. 



Stiffness and damping coefficient matrices, 
where i,j= x, y, a, /9. 


k ii - K, (C/P 0 R 2 ) 


dimensionless stiffness coefficient matrix, where 
i,j= x, y, a, p. 


L 


seal length. [L] 


M*. M y 


components of fluid film force about x and y 
axes. [F-L] 


m = M/(P 0 R 3 ) 


dimensionless fluid film moment. 



h 

P 

P = P/P 0 

P, P r 
P p . P s 
Po 

Q r 

q r = Q r (12jx/P 0 C 3 ) 

R 

Re* = ph^p/zi 2 
Re 0 ‘ = pC 3 P 0 /(Rzx 2 ) 
t 

u = U (13 xR/C 2 P 0 ) 

v = V (1 ?ziR/C 2 P 0 ) 
u, V 

Ub-Uj 

X.Y.Z 
z = Z/R 
a = A (2L/C) 


unit vector normal to fluid film boundary. 

local pressure. [F/L 2 ] 

dimensionless local pressure. 

Left and right boundary pressures. [F/L 2 ] 

Pocket and supply pressures. [F/L 2 ] 

Reference pressure, used for scaling the 
pressure field, which is normally set equal to P s , 
Pp. P, or P r . [F/L 2 ] 

flow from pocket or recess. [L 3 /T] 

dimensionless flow from pocket or recess. 

seal radius. [L] 

local Reynolds number based on pressure- 
driven flow. 

refernce Reynolds number based on pressure- 
driven flow. 

time. [T] 

dimensionless circumferential component of fluid 
velocity. 

dimensionless axial component of fluid velocity. 

circumferential and axial fluid velocity 
components, averaged across the film. [L/T] 

linear velocity of housing and journal surfaces 
(equal to R&> b , Rwj, respectively). [L/T] 

cartesian coordinates. [L] 

dimensionless axial coordiante. 

misalignment ratio about the x-axis. 
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0 = B (2L/C) 

6 x- € y 

e = e/C 

e 

A b = euU„R/(C 2 P 0 ) 

A, = 6f I U J R/(C J P 0 ) 

A ,= P C 6 P c / {288 A..VP 2 ) 
A. - K, (Re 0 'C/288 R) 

M 

P 

“b- 

r = t (C 2 P 0 /13xR 2 ) 


misalignment ratio about the y-axis. 
components of rotor eccentricity ratio, 
rotor eccentricity ratio, 
circumferential coordinate, [radians] 
dimensionless velocity of housing surface, 
dimensionless velocity of rotor surface, 
coefficient of orifice restriction, 
coefficient of pressure drop at inlet to film, 
fluid dynamic viscosity. [F-S/L 2 ] 
fluid density. [F-T 2 L' 4 ] 

angular velocity of housing and journal surfaces. 
[rad/T] 

dimensionless time. 


4-60 



5.0 INDUSTRIAL COMPRESSIBLE CYLINDRICAL CODE, GCYL 


The program GCYL (acronym for gas lubricated cylindrical) is used for analyzing 
a variety of seals that can be defined in a cylindrical coordinate reference 
frame. Figure 5-1 shows solid ring configurations and Figure 5-2 shows typical 
sectored ring configurations that the program analyzes. The capabilities of the 
program include the following: 

e Varying geometries, as indicated on Figures 5-1 and 5-2. 

e Variable or constant grid representation. Maximum grid size is 30 
grid points in the axial direction and 74 grid points in the 
circumferential direction. Figure 5-3 shows a typical grid network. 

The circumferential parameter is 0, and the axial parameter is Z. 

The grid points are identified in the axial direction as I and in 
the circumferential direction as J. The extent of I is 1-*M, and the 
extent of J is 1-*N. 

e Specified boundary pressures or periodic boundary conditions in 
the circumferential direction. 

e Axial symmetry option 

e Four degrees of freedom, x and y translations of rotor origin and 
angular displacements about the x and y axes through the rotor 
origin 

e Determining load as a function of shaft position or determining 
shaft position to satisfy a given load. 

e External Pressurization (Hydrostatic) of inherently compensated 
orifices, spot recesses or full recesses. 

e Choice of English or SZ units 

The output of the program includes: 

e Clearance distribution 

e Pressure distribution 

e Leakage along specified flow paths 
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Figure 5-1 Leakage Path Geometries (Floating Ring) 











• Load and Load angle 


• Righting Moments 

• Viscous dissipation 

• Cross -coupled, frequency -dependent, stiffness and damping 
coefficients 

• Plotting routines for the pressure and clearance distribution 

The program has been written for a PC environment using OS/2 as an operating 
system. Relatively large dimensions have been utilized which would exceed the 
memory limitations of a DOS environment. The FORTRAN code however, would be 
amenable to other systems that employ FORTRAN 77 as long as memory is sufficient. 
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5.1 THEORETICAL DESCRIPTION AND NUMERICAL METHODS 


General Theory 


Reynolds' equation for compressible flow for a cylindrical geometry* is as 
follows: 


_L JL 

r2 39 




+ 12y^ 


(5-1) 


The equation is made dimensionless with the following definitions. (Upper case 
variables are dimensionless). 

Z - z/R, H - h/c 0 , T - t/t 0 , P - p/p o . 


6vwR2 

m _ t m 


p c ^ 

*o o 


12uR z 

p C 2 
r o o 


Substituting the dimensionless variables into Reynolds' equation produces a 
dimensionless Reynolds' equation. 



(5-2) 


For steady-state solutions, the time dependent term on the right hand side is 
eliminated except for the computation of spring and damping coefficients. 

In the solution methods subsequently described, the Reynolds' equation is not 
applied directly. The Reynolds' equation represents the divergence of the mass 
flow at any grid point. The more convenient cell method is to conduct a mass 
balance directly, and not the divergence of the mass flow at each point. 

Formation of Equations for Determining Pressure Distribution 


*Nomenclature is included in Section 5-5 





The general method of solving for the pressure distribution is the cell 
method^ whereby a flow balance through a cell volume is accomplished. The 
perimeter of the cell extends halfway between the grid point and its four neigh- 
boring points. A typical cell is shown by the dashed lines on Figure 5-4. The 
principal grid point is at Row i (length direction) and Column j (circumferen- 
tial direction). For convenience of programming the grid points are numbered 
for each cell sequentially from 1 to 9 with grid point 5 being the principal 
point. The corners of the cell boundaries are also numbered from 1 to 4. 

Figure 5-5 shows the flow balance through the cell. There are eight flows 
across the cell boundaries, and there can also be a source (or sink) flow into or 
out of the cell control volume. The reason eight flows are used in lieu of four 
is that it permits discontinuous clearance boundaries at grid lines (such as 
Rayleigh steps) without taking derivatives across a discontinuous boundary. 


The net flow through a cell can be expressed as: 

+ 


<12 


<34 


AZ i . n - 

AZ i-l 

2 + Q 12 

2 

AZ. 

AZ, , 

i - 

~ 2 ~ ~ Q 34 

i-1 

2 


+ Q 


- Q 


AO 

14 2 

AO 
2 


AO 


23 


+ Q 
- Q 


i-i 


14 2 

AO 
23 


(5-3) 


izl - 


‘in 


Ql2 + means the mass flow per unit length across the plus side of cell boundary 
1-2, etc. 

The Q's are dimensionless mass flows per unit length, except for Q{ n which is a 
dimensionless source inlet flow. 


In the 6 direction 


Q 


If 



paf 


In the length or Z direction 

, 3P M 
Q - -PH 3Z 2 

♦♦References are identified in Section 5-4 
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where Q is defined as 


Q 


12y G c T a 

p 2 c 0 3 

r o o 


( 5 - 6 ) 


(Primed values of P are absolute pressures; unprimed values are gage pres- 
sures). 


An optional flow can enter the cell from an external source, which can be treat- 
ed as an inherently compensated orifice, or a conventional orifice restriction. 
Inherent compensation presumes the orifice area is the surface area of a cylin- 
der circumvented by the hole size and length equal to the clearance under the 
inlet hole. The conventional orifice area is the area of the hole. The conven- 
tional orifice generally discharges into a recess that allows the flow velocity 
to dissipate into a region of constant pressure. Two types of recesses are 
permitted; a spot recess which is treated as a source at one grid point, or a 
recess of finite length in the axial and circumferential directions, which is 
fed by an inlet orifice. 


Pressures are taken as the average pressure across the boundary. For example: 

P i’j + P iJ +1 


P 12" 


and 


( 5 - 7 ) 


3P 

36 


p i,j+l “ P i»J 


12 


66 , 


etc. 


( 5 - 8 ) 
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By substituting the pressures and pressure derivatives (equations 5-7 and 5-8) 
into the mass flow balance equations (equation 5-3 and 5-4), an equation is 
derived that is a function of the five pressures, ?2f ^ 4 > ^ 5 1 an< i ^8* anc * t ^ ie 
clearances taken at the cell corner points Hj — H 4 . Each cell corner point film 
thickness is computed in* the clearance routine by appropriate values of Z and 9 
and is designated as as HC{, i = 1,4. For example, HC^ is the clearance at the 
cell corner point 1 . 


An optional flow can enter the cell from an external source, which is treated as 
an inherently compensated orifice or the usual hole size orifice restriction. 
Point sources pose numerical instability problems , which are circumvented by 
applying fine grids surrounding the source points. Flow through the orifice is 
given as: . . - .^ 1/2 


<in 


OFC x AO x 


ir h 


:ni 


(5-9) 


where, 


12 V C d / 2YG c T a 
OFC - _ 7-7 . / 


pTT 

*0 o 


(5-10) 


! irdoH 5 for Inherent Compensation 

Trdo^ for orifice compensation 
4 (spot recess or full recess) 


(5-11) 
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where CR 


f 2.0 


[Ty+I)J 


(£) 


_J —1 

Also, if F R >1.0, R 

71 71 




end P S 


Pi 


R 


(5-13) 


(5-14) 
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This condition implies backflow through the orifice. 


The primed values indicate absolute pressures (i.e., P'r = Pr+1). 

Thus, a system of numerical equations can be derived as a function of five pres- 
sures. There is an equation for every grid point. 

f (P r P 2 P 5 > i,j - 0 (5-15) 


The system is non-linear since it is dependent upon multiples of P and its 
derivatives. 


The solution process starts by assuming a pressure distribution, and using 
Newton-Raphson iteration until the functions f converge to zero within a 
prespecified truncation error. In equation form, the iteration process is: 


i.J 


(old) + Z 


k-1 


nJ° ld) 

,P K 


p (new) _ p (old)j 

K K i 

I 


(5-16) 


where the partial derivatives are explicitly determined, e.g. 


3P 
0 K 


f(P^» ^ 2 * •••^5 i 

• f (P^ t ^2* * * * * c/2» ••• ^ j 


(5-17) 


The actual convergence is not on f, but on p^^ new ^ - for when the 

difference vanishes* the condition that f = 0 is satisfied* 
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The Column Method Solution of Newton-Raphson Equations 


The column method*- J is used to solve the new pressures in the set of mxn 
equations defined by equation 5-16. The advantage of the column matrix method 
is that its inversions are M x M rather than M x N so that its use saves computa- 
tional time. 

The linearized N-R equations may be written in the form: 

c j p j + E j p j-l + D j p j+l * R j (5-18) 

For each value of j, Pj is a vector containing the jth column of new pressures , 
Rj is the right hand side column vector and Cj, Ej and Dj are in general 
tri-diagonal matrices. 

Case 1 - Pressure Prescribed at Start and End of Pads 


Equations of form 5-18 are written at all points in the grid corresponding to i = 
1, 2, +M and j * 2, 3, **N-1 with boundary column vectors and prescribed. 
Look for a solution in the form: 

p j_l = Aj Pj + Bj (5-19) 

Where Aj is an M x M matrix and Bj is a vector. Use equation 5-19 to eliminate 
Pj-l appearing in equation 5-18. 

(Cj ♦ Ej Aj) Pj + Ej Bj + Dj Pj +1 = Rj (5-20) 

Solve equation 5-20 for Pj to obtain: 

Pj = “I j Dj Pj+i+ I j(Rj — Ej Bj) (5—21) 
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(M x M matrix) 


Where Ij = (Cj + Ej Aj ) -1 
Set j = j+1 in equation 5-19 to obtain 

P j * A j + 1 P j + 1 + B j+1 (5-22) 

Compare coefficients in equation 5-21 and 5-22. 

Aj +1 = — Ij Pj» Bj+i = Ij (Rj — Ej Bj ) (5—23) 

Set A 2 *0, B 2 * Pi 

Use equation 5-23 to compute A 3 , A 4 , — , A^ and B 3 , B 4 — Bn. 

Since Pfl is given and all Aj and Bj are computed, we may use equation 5-19 to 
compute Pn-1» Pn-2> p N-3 » > p 2 * 

Review of General Procedure for Non-Periodic Boundaries 

1) Set A 2 = 0 

B 2 = Pi 

2) Compute Aj+i, Bj+i 

Aj +1 * "Ij Dj 

B j+i = Ij (Rj - Ej Bj) j - 2, N-l 

where Ij * (Cj + Ej Aj)”* 

3) Compute Pj 

Pj_l = Aj Pj + Bj j * M, 2 

Case 2 - Column Method for Periodic Boundaries 

Pj, Bj, Rj, Zj are vectors. N' = N -1 
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For periodic boundaries, the condition is that Pi - Pn* At the boundary, j“l, 
the general equation is: 


Ci Pi + Ei P N ' + Di P 2 = Rl 


( 5 - 24 ) 


At column N', the equation becomes 

Cm' Pn’ + Eh* Pn'-i + Dn’ pi ■ Rn' 

To satisfy the boundary conditions, a solution is assumed of the form: 


( 5 - 25 ) 


Pj-1 - Aj Pj + Bj + Fj P N ’ 


Ai * 0, Bi = 0, Fi * 6 (Kronocker delta matrix) 


( 5 - 26 ) 

( 5 - 27 ) 


Returning to the general equation: 


Cj Pj + Ej Pj— 1 + Dj Pj+1 Rj 

Substituting for Pj-i from equation 5-26, the following results: 


( 5 - 28 ) 


(Cj + Ej Aj) Pj + Ej Bj + Ej Fj P N ' + Dj Pj+i = Rj 


Ij = (Cj + Ej Aj) 


-1 


( 5 - 29 ) 

( 5 - 30 ) 


Then, 


Pj = “ Ij Dj Pj+1 + Ij (Rj “ Ej Bj) — Ij Ej Fj Pn* 


( 5 - 31 ) 


From equation 5-26: 


Pj “ Aj+i Pj+l + Bj+i + Fj+i P N ' 


( 5 - 32 ) 


Comparing equations 5-31 and 5-32: 


Aj+i = “Ij Bj, Bj+i = Ij (Rj “ Ej Bj), Pj+1 — “Ij Ej Fj j— 1» 2, — , N— 1 


(! 


For Pn = Pl> we obtain from equation 5-26: 


- 33 ) 
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P N ' = An'+i Pi + Bjj'+i + Fn' + 1 Pn’ 
After rearranging: 

p N « = (6 - f n »+i) _1 (An'+i p i + b n'+i> 


(5-34) 


(5-35) 


or PN' = y N' p l + Z N' 


(5-36) 


where 


Yn 1 * (« - FwUir 1 A N * + 1, Z N f * (5 - Fn' + iT 1 B N *+i (5-37) 


Substituting equation 5—36 into 5—26 we obtains 


p N'-l * A N ' (Y N * Pi + Z N ’) + B N * ♦ Fjj' (Y N . Pi + z N .) 


= (A N « Y N ' + F N » Y N ') Pi + a n * z n * + b n « ♦ F N ' Z N ' (5-38) 


= Yy'-i pi + zn'-i 


where 


Yn'-1 = An' Yn' ♦ Fu» Yn', Zn'-1 = An' Zn' + Bn' + Fn’ Zn' (5-39) 


Similarly, 


Pn'-2 s Ajj'-I (Yn' - 1 Pl + Z N '-l) + Bn'-1 + Fn'-1 (Yn' Pi ♦ Zn') (5-40) 
= (An'-i y N '-i + f n'-i y n') pi + a n*-i z N '-i + b N '-i + f n *-i z n * 


= Yn' - 2 p l + z N'-2 


(5-41) 


Yj-1 = Aj Yj ♦ Fj Yn' 

Zj-1 = Aj Zj ♦ Bj + Fj Z N ' 


(5-42) 
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Therefore, in general: 


P j -1 = Yj_i Pi + Zj-i or Pj = Yj Pi + Zj (5-43) 

Pi = (5 - Yi) _1 Zi (5-44) 

Review of General Procedure for Joined or Periodic Boundaries 


1) 

Compute Aj+i 

» Bj+l, Fj+i 



Aj+l s 

‘ -Ij Dj 



Bj+1 s 

■ Ij (Rj -Ej Bj) 

j-1. N-l 


Fj + 1 - 

= -Ij Ej Fj 




Ai =0 

1 j » <Cj + Ej A 



Bi =0 




F i =6 


2) 

Compute Y ^' f 

Zn' 



y n * = 

(6-F N ) _1 A N 



Z N * = 

(6-F N )’ 1 B N 


3 ) 

Compute 




Yj-l - 

! Aj Yj ♦ Fj Yjj' 

j = N' -*• 2 


Zj-i - 

! A j Zj + Bj + Fj 

z N ' 

4 ) 

Compute ?i m 

(6 - Yi)" 1 Zi 


5 ) 

Compute Pj « 

Yj Pi ♦ Zj 

j = 2 ♦ N' 
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The coefficient matrices Cj, Ej and Dj, and the right hand side vector Rj, are 
easily formulated. Cj contains all the coefficients multiplied by Pj. By exam- 
ining equation 5-16, it is seen that for any row i and column j that values of C 
are: 


c i, 

i-l» j = 

3£ 5 

3*8 

r • 

II 

**"» 

•H 

3*5 


3*5 

C i, 

i+lt J " 

3*5 


(5-45) 


Similarly the coefficient matrix Ej contains the elements: 


and 


J 


D i, i. J 


3 f 5 

3*4 

111 

3*6 


(5-46) 

(5-47) 


Rj contains all elements not multiplied by the pressure 

5 

Rj - -fij ( ° ld) + Z till P K (° ld) 


1 3P, 


K 


(5-48) 


Separate subroutines are used depending upon the pressure boundary conditions. 
The subroutine COLP implements the column method for prescribed boundary condi- 
tions while COLJ does it for periodic or joined boundaries. The subroutine 
COEFC forms the C matrix coefficients while COEFF forms the coefficient matrices 
and right hand side vector D, E, R respectively. 

Film Thickness Distribution (see Figure 5-6) Eccentricity and Misalignment 


In vector format, the clearance due to eccentricity and misalignment at any 
angle 0 and at distance z' from the mid-plane is: 


h ■ 

[C o e r - e x i - e y j - oi x z'k - BJ x 

»*) s 

(5-49) 

h - 

C - e cosfl - e sin6 + az' sin0 - 

A V V 

8 z’cosO 



ox y 


(5-50) 

- 

C - (e + 0z') COS0 - (e - az') 
o ' x y 

sin0 
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Using dimensionless variables, equation (5 - 50) becomes! 
H - 1 - |e x + 6 (Z - L/2)R j cos6 

C o 

- |e - a (Z-L/2)R \ sin6 


which is set equal to: 


H = 1- (e x + eg) cos 0 - (e y ♦ e 0 ) sin 0 


where 

e - 8 (Z ~ L/2)R 


- « (Z - L/2)R 
° C o 


(5-51) 


(5-52) 


(5-53) 


Preloaded Seals 


Preloaded seals (see Figure 5-7) can be modeled by adding an additional eccen' 
tricity in the x and y directions. 


e 


e 


X 

■ COS0 

PR 

PR p 

y 

« E__ sin6 

PR 

PR p 


(5-54) 


where 


e p * * x eccentricity due to preload 
e p y ” y eccentricity due to preload 
8p ■ preload angle 


Rayleigh Step 


The grid network for the Rayleigh step is shown on Figure 5-8. The boundaries of 
the step are defined by the lower left and upper right corners of the depressed 
region. Interior grid points include the step height in the clearance distrib- 
ution. . 
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X 


Keyword 

Variable 

Description 

START 

THETAS 

Pad Start Angle 

PADANGLE 

THETAP 

Pad Angle 

PIVOT 

PVT 

Pivot Angle 

PRELOAD 

EPR 

Offset/Clearance 


Figure 5-7 Preloaded Seal 
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Figure 5-8 


Raylelgh-S.tep 
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Axial Taper 


An axial taper is indicated as Figure 5-9. If Z > Z t then 


H = H + 6(Z-Z t ) 


(5-55) 


5-22 




Power Loss (Torque) 


Power loss is obtained by integrating the viscous shear forces across the film. 
From Figure 5-10, a force balance on an element produces: 


12.3T (5-56) 

3x 3z * 


but 



(5-57) 


Therefore, 


»E 

3x 


3 ? 


(5-58) 


Integrating, 


3U 

3z 


1 l£ 
U 3x 


z + C 


1 


U 


i iE 
U 3x 


z? 

2 


+ C^z + C 


2 


(5-59) 
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The boundary conditions are: 

U - 0 z - 0 ••• c 2 ■ 0 


U - U when z ■ h 


Substituting: 


U - - Jl 2 + c ,h 


W 3x 2 


(5-61) 


(5-62) 
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Figure 5-10 Viscous Power Loss 
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Therefore, 




(. 5 - 63 ) 

( 5 - 6*0 

( 5 - 65 ) 

( 5 - 66 ) 

( 5 - 67 ) 

( 5 - 68 ) 

( 5 - 69 ) 

( 5 - 70 ) 
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Computation of Seal Flows 


The program computes the flow across specified axial and circumferential grid 
lines. A total of four grid lines can be prespecified. The subroutine FLOCIR 
determines flow across a circumferential line and the subroutine FLOAXL computes 
the flow across an axial grid line. 

Circumferential Flow Line (see Figure 5-11) 


There are three types of points to consider. A point on a grid boundary J -1 or 
J =N, and an interior point. Also, a flow line on I =H requires special consid- 
eration. For each point, a flow balance is accomplished through the cell 
surrounding the point as depicted on Figure 5-11. Consider an interior grid 
point on an interior grid line (I = M). 


where 



DZj/Z 


(5-71) 

(5-72) 

(5-73) 

(5-74) 


The remaining flow components are similarly computed and Qc(I,J) determined. 

1L[ . , III 

At J * 1, 30 I is computed by forward difference and is equal to 3 q J • 

34 + + 12 

The pressure P 34 * Pl2 • 

The clearance H 4 is not a regular grid point clearance and thus is not included 
in the grid clearance array. H 4 is computed as the average of Hj.j and Hi+i, j. 

The grid line mass flows are accumulated to obtain the total flow across the 
grid line. 
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Similar procedures are employed for computing flows across axial lines (see 
Figure 5-12) • 
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FREQUENCY DEPENDENT SPRING AND DAMPIN G COEFFICIENTS 


Discretization has been carried out with the use of the cell method [1] , which 
involves a flow balance about each grid point. 


fV-QdA - jQ-ndS « (1 +P) HdA 


(5-75) 


where Q - the mass flow vector per unit length. 

The equality of the first two terms comes from the divergence theorem. 
In numerical format the right hand side becomes 




(5-76) 


where, 

Alj - 1 (AOj+Ae^) (AZ^AZ^) (5-77) 

Generally, a small perturbation analysis is used for determining frequency 
dependent spring and damping coefficients and solving the complete equation (5- 
75). A small perturbation analysis, however, is generally limited to concentric 
operation and produces complex expressions for the perturbation coefficients. 
Identical results can be achieved by direct numerical perturbation of the 
difference equations used in the column matrix solution approach. This method, 
which is described below, avoids algebraic error in determining the perturbation 
coefficients and may be used in complex situations where analytical determination 
of the perturbation coefficients is not feasible. 


After desired convergence of the Newton-Raphson process has been achieved under 
steady (unperturbed) conditions, the resulting steady state pressure vectors are 
denoted as (P) and the coefficient matrices as [C 3 ], etc. and as before the 
steady state equation becomes: 


The eccentricity components can be perturbed individually by an amount r), and the 
matrix [C J ] recalculated at the new film thickness (but old pressure 
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distribution, P) ; then subtract [C J ] at the old film thickness and divide the 
difference by T] to numerically obtain the partial derivative of [C^] with respect 
to the eccentricity perturbation. This partial derivative will be denoted by 
[C J,k ] . Thus, 




[gju - ML, 

n 


(5-79) 


The matrices [E J>k ], [6^] and {R J,k ) are obtained in a similar manner. Equation 
(5-78) may now formally be differentiated with respect to € k to obtain the 
expression: 

[<?W) *WA) 

(5-80) 

where {P’ k j) - 3{Pj)/d€ k is the zero frequency stiffness pressure. This expression 
does not yet contain the time dependent terms found on the right hand side of 
equation (5-75). It is assumed that a sinusoidal disturbance is applied to the 
shaft, such that the clearance and pressure derivatives are affected as follows: 

H * e io1 ; - Pj k e ioB (5-81) 


To complete the process the right hand side of equation (5-75) is differentiated 
with respect to e k and the results added to the right hand side of equation (5- 
80) with d/dt replaced by i a. The terms to be added to the right hand side of 
equation (5-80) in this manner are -ia[0*] {Pj’ k } -i<7{R J,k ) where [C J ]are diagonal 
matrices whose components are 

C ii * (5-82) 

Because a cell can have clearance discontinuities, such as a step, it is 
advantageous to partition the cell into 4 components as indicated on Figure 5-2, 
and then sum the components to obtain [£■*]. Thus equation (5-82) becomes: 

c j ii = HC 1 A 1 + HC 2 A2 +HC 3 A 2 + HC 4 A i (5-83) 


where; HC X is the clearance at the comer point 1 of the cell and 
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*1 


4 


a 2 


(Ae^AZ,.,) 

4 


; etc. 


( 5 - 84 ) 


and (R J,k ) are column vectors whose components are 


dH ±i 


Xi'* * 


( 5 * 05 ) 


By combining terms, the final set of linear difference equations for the complex 
stiffness pressure derivatives {P' k j) are obtained 

[c-J]{Pi*} ♦[**]{*&} +[6%p&) - {« J *} -[^■TO -\F" V,.i} 


where , 


[C*J] = [^l + iotc^ ; = {£*'*} 


(5-86) 

( 5 - 87 ) 


The system of equations given by equation (5-86) is solved by the column method 
in a directly analogous manner to that used in solving the steady state equation. 
The principal difference is that all the matrix operations are performed using 
complex arithmetic. Integration of the real part of the pressure derivatives 
yields the stiffness while the complex parts when integrated and divided by a 
yields the damping. 
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5.2 SAMPLE PROBLEMS 


5.2.1 Ravleieh- Step Seal 

The first sample problem is a four pad Rayleigh-step seal (refer to Figure 5-8). 
The geometry and operating conditions are as follows: 

• Number of pads - 4 

• OPTION - 1, which means the position of the seal will be pre-specified. 

• Seal length - 3.852 in. and the symmetry option will be used. 

• Variable grid will be used in the axial and circumferential directions . 
Since symmetry has been applied in the axial direction, the variable grid 
length equals half the actual length, and is equal to 1.926. in. 

• The grid will be made finer at the step boundaries where sharp pressure 
gradients are expected to occur. 

• Seal diameter — 1.9685 in. 

• The step height is 0.00165 in. deep and is located at the leading edge 
of the pad, 5 degrees from the x axis and the lower left corner of the 
step is 0.655 in from the inside radius. The end of the step is 70.6 
degrees from the x axis, and since symmetry has been invoked the axial end 
of the step as represented on the grid is 1.926 in from the inlet end, or 
at the end of the grid. 

• The specific heat ratio of the lubricant is 1.4. 

• The gas constant is 250,000 in 2 /(s 2 -°R) 

• The absolute temperature is 530°R 

• The absolute viscosity is 3.0 x 10‘* lb-s/in 2 

• The eccentricity ratio - 0.2 

• The eccentricity angle is 270° 

• The shaft speed is 70,000 rpm 
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• The reference pressure is 200 psig 


• The boundary pressures are all 0 psig, or 200 psia. 

Results of the problem are indicated on Table 5-1. 

Figures 5-13 and 5-14 show the clearance distribution and the pressure 
distribution produced by the plotting programs. These plots clearly show the 
highly loaded pad, which is pad number 3 (highest pressure level and lowest film 
thickness level) . 
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TABLE 5-1 


Summary of Results 
Sample Problem 1 


-JOURNAL l LOAD POSITION 
ECCENTRICITY 
ECCENTRICITY ANGLE 
MINIMUM FILM 
LOAD 

LOAD ANGLE 


.40000 
-90.00 DEG 
.0006015 IN 
27.54 LB 
61.44 DEG 


POWER LOSS 


.4555 HP 


LEAKAGE AT 1 ■ 1 * -.44783E-03 LB/S 
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Figure 5-13 Rayleigh Step Seal - Clearance Distribution 
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PRESSURE DISTRIBUTION 


3CYL 

SAMPLE CASE 1 : RAYLEIGH STEP SEAL 

DI !™ETER = 1.969 IN SPEED = 70000.00 RPM 

LENGTH = 3-852 IN 

CLEARANCE = .001000 IN 



Figure 5-14 Rayleigh Step Seal Pressure Distribution 


174E+02 
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5.2.2 Sample Problem 2 - Non-grooved Lobe Seal 

The non-grooved lobe seal is characterized by offset lobes that are joined at 
their apexes in a continuous fashion as opposed to a lobe seal where the lobes 
are separated by axial grooves. Such a seal is depicted on Figure 5-15; it would 
be manufactured by a broaching process . To analyze this type of seal with the 
GCYL code the key word SECTOR must be invoked followed by the number of sectors, 
the lobe preload and preload position within the lobe (see Figure 5-7 for 
definition of preload). For this example, a lobe hydrodynamic geometry was 
combined with a external pressurization through source points at the midplane of 
the seal. The geometry and operating conditions are as follows: 

e Seal Diameter - 2.25 in. 

e Seal Length - 1.625 in. 

• Seal reference clearance - 0.0005 in. The reference clearance is the 
clearance prior to preload. 

• Number of pads - 1. A sectored seal is always considered as a continuous 
seal although discontinuities exist in the clearance distribution. Thus, 
the n umb er of pads are always unity and the JOINED option is always 
applied. 

e The preload on each lobe is 0.5 which means at the pivot position the 
lobe is eccentric toward the shaft one half of the reference clearance 
(see Figure 5-7). 

e The pivot angle of the first sector is 150° from the x axis, and since 
the first lobe is 90° from the x axis the pivot position is located at the 
mid angle of each lobe. 

• The viscosity of the gas is 3 x 10' 9 lb-s/in 2 

• The gas constant is 2.5 x 10 5 in 2 /(s 2 -°R) 

• The ambient temperature is 510°R 

• The total number of orifices is 27, 9 in each sector, located at the 
midplane of each sector. One orifice is located at each interior grid 
point at the mid plane of the bearing. 
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• The orifice diameter is 0.015 in. 


• The coefficient of discharge of each orifice - 0.9 

• The supply pressure to the source orifices is 120 psig 

• The operating speed - 70,000 rpm 

• The reference pressure is 14.7 psig 

• The pressure along the boundaries is 0 psig 

Table 5-2 summarizes output data. Since the shaft is concentric within the seal, 
total load capacity is zero. The most important information is the leakage flow. 

The clearance and pressure distribution are shown on Figures 5-16 and 5-17 
respectively. Notice the discontinuities in the clearance distribution because 
of the lobed geometry. The proximity of the source points to each other makes the 
pressure distribution appear as a line source. 
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TABLE 5-2 


Summary of Performance 
Sample Problem 2 


'JOURNAL t LOAD POSITION 


ECCENTRICITY 

* 

.00000 

ECCENTRICITY ANGLE 

■ 

.00 DEG 

MINIMUM FILM 

* 

.0002500 IN 

LOAD 

m 

.9904E-12 LB 

LOAD ANGLE 

m 

•53.73 DEG 

POWER LOSS 

m 

1.233 NP 

LEAKAGE AT I » 1 

m 

-.U924E-03 LB/S 

LEAKAGE AT I - M 

m 

.H924E-03 LB/S 

RIGHTING MOMENT 



ABOUT X-X 

HX ■ 

-.2814E-13 LB-IN 

ABOUT Y'Y 

HY ■ 

.1784E-14 LB-IN 
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FILM THICKNESS DISTRIBUTION 


GCYL 

SAMPLE CASE 2: SECTORED SEAL 

DIAMETER = 2.250 IN SPEED = 70000.00 RPM 

LENGTH = 1 .625 IN 



Figure 5-16 Clearance Distribution, Sectored Lobe Seal 


5-43 





Figure 5-17 Pressure Distribution, Sectored-Lobe Seal 
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5.2.3 Sample Problem 3 - Three Lobe Seal 


This problem deals with the hydrodynamic portion of a 3 -Lobe seal where the lobes 
are separated by axial grooves. Figure 5-7 shows the general geometry and key 
parameters. The principal parameters are the preload and pivot angle. The 
following are geometry and operating conditions: 

• OPTION -2 which means the position of the seal to satisfy a given load 
will be determined. 

• International units apply; parameter SI invoked. 

• Stiffness and damping are to be calculated in two degrees of freedom, x 
and y, at an imposed frequency equal to running speed of 50,000 rpm. 

• The number of pads -3 

• The start of the first pad is at 100°, and the pad extent is 100° 

• The pad preload is 50 Z of the reference clearance, and the preload for 
the first pad occurs 150° from the x- axis, which means the preload is in 
the center of the pad. 

• The shaft diameter is 0.0508 m. 

e The hydrodynamic length is .0254 m 

• The reference clearance is 1.27 x 10' 3 m. 

• The lubricant viscosity is 2.07 x 10* 3 N-s/m 2 
e The absolute temperature is 283°K 

e The ratio of specific heats of the gas is 1.4 
e The gas constant is 290.32 m 2 /(s 2 -°R) 
e Symmetry is applied in the axial direction 
e The load to be supported is 200.16 N 

e The angle at which the load is applied is 270° from the x-axis. 
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• The initial eccentricity guess is 0.2, and the initial guess on the 
displacement angle is 270° from the x-axis. 

• The shaft speed is 50,000 rpm 

• The reference pressure is 8.274 x 10 5 Pa. 

• The boundary pressures are all 0 gage. 

Table 5-3 indicates the steady-state performance and stiffness and damping 
coefficients. 

Graphical representations of clearance and pressure distributions are shown on 
Figures 5-18 and 5-19 respectively. 
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TABLE 5-3 


Summary of Performance 
Sample Problem 3 


-JOURNAL l LOAD POSITION 
ECCENTRICITY 
ECCENTRICITY ANGLE 
MINIMUM FILM 
LOAD 

LOAD ANGLE 


POWER LOSS 

■ 

186.1 

If 

LEAKAGE AT I • 1 

■ 

- • 12875E-04 KG/S 

STIFFNESS COEFFICIENTS 



PRINCIPAL X 

KXX « 

.1460E+09 N/M 

CROSS-COUPLED 

KXY ■ 

- .23796+08 N/M 

CROSS-COUPLED 

XXA « 

.0000 

N/RAD 

CROSS-COUPLED 

KXB » 

.0000 

N/RAD 

CROSS-COUPLED 

KYX ■ 

•.309OE+OS N/M 

PRINCIPAL Y 

KYY • 

.10026+09 N/M 

CROSS-COUPLED 

KYA - 

.0000 

N/RAD 

CROSS-COUPLED 

KYB « 

.0000 

N/RAD 

CROSS-COUPLED 

KAX « 

.0000 

N-M/M 

CROSS-COUPLED 

KAY « 

.0000 

N-M/M 

PRINCIPAL A 

KAA ■ 

.0000 

N-M/RAD 

CROSS-COUPLED 

KAB > 

.0000 

N -M/RAD 

CROSS- COUPLED 

KBX - 

.0000 

N-M/M 

CROSS-COUPLED 

KBY ■ 

.0000 

N-M/M 

CROSS-COUPLED 

KBA ■ 

.0000 

N-M/RAD 

PRINCIPAL B 

KBB ■ 

.0000 

N-M/RAD 

DAMPING COEFFICIENTS 



PRINCIPAL X 

DXX ■ 

9939. 

N-S/M 


CROSS- COUPLED 

CROSS-COUPLED 

CROSS-COUPLED 

CROSS-COUPLED 

PRINCIPAL Y 

CROSS-COUPLED 

CROSS-COUPLED 

CROSS-COUPLED 

CROSS-COUPLED 

PRINCIPAL A 

CROSS-COUPLED 

CROSS-COUPLED 

CROSS-COUPLED 

CROSS-COUPLED 

PRINCIPAL • 


DXY 

DXA 

DXB 

DYX 

DYY 

DYA 

DYI 

DAX 

DAY 

DAA 

DAS 

DBX 

DBY 

DBA 

DBB 


•22103 
129.38 DEG 
.0000037 M 
200.2 N 
•90.00 DEG 


-.1089E+05 N-S/M 
•0000 N-S/M 
•0000 N-S/RAD 
5140. N-S/M 
•1393E+05 N-S/M 


.0000 

.0000 

.0000 

.0000 

.0000 

•0000 

.0000 

.0000 

.0000 

.0000 


N-S/RAD 

N-S/RAD 

N-M-S/M 

N-M-S/M 

N-M-S/RAD 

N-M-S/RAD 

N-M-S/M 

N-M-S/M 

N-M-S/RAD 

N-M-S/RAD 
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Figure 5-18 Clearance Distribution — Three— Lobe Seal 
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Figure 5-19 Pressure Distribution - Three-Lobe Gas Seal 


5-49 



5.2.4 Sample Problem, 4 - T-Shaped Sectored Seal 


This problem deals with an actual helium buffered seal analysis and design (for 
SSME) that was accomplished for NASA. A design that incorporates a self - 
adjusting clearance that can accommodate thermal and centrifugal distortions and 
shaft dynamic excursions avoids many of the problems associated with captured 
clearance designs. The sectored ring seal provides the desired self adjusting 
clearance features. The general configuration of the sectored seal is shown on 
Figure 5-20. The sectors consist of T-shaped sections mated to each other at each 
end with sealed joints. The sectors can move relative to each other 
circumferentially and that is how the seal accommodates variations in the sleeve 
dimensions due to thermal expansions and contractions and centrifugal growths. 
The T-shaped sector was chosen because it is a symmetrical shape and the various 
fluid and friction forces can be designed to avoid upsetting moments on the 
individual sectors. An overlapping V joint prevents a direct clearance path 
between the hydrogen and oxygen ends of the seal. Each sector is supported by a 
hydrostatic fluid- film on its inner circumference and along the side walls which 
forms a friction free secondary seal to permit free radial motion of the sectors 
in response to sleeve movements. The fluid-films are predominantly hydrostatic 
to avoid any pitching tendencies introduced by the hydrodynamic effects. The 
hydrostatic bearings are fed by the buffer pressure on the outside diameter of 
the seal. Figure 5-21 shows the pressure distribution and force balance on the 
individual sectors. 

This sample problem describes one case conducted in the analysis of the 
circumferential hydrostatic seal on one of the sectors. The geometry and 
operating conditions are as follows: 

• The number of pads to be analyzed is 1. 

• OPTION - 2, which means the position of the sector to satisfy a given 

load will be determined. 

• The load applied is 370 lbs. 

• The load angle from the x -axis is 270° 

• The initial guess on the eccentricity of the seal is 0.5 

• The initial guess on the eccentricity angle is 90° 

• Variable grids are used in both the axial and circumferential direction. 

The grid is made very fine around the source points. The starting angle of 
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the sector is 30° from the x-axis and its angular extent is 120°. the axial 
length of the seal is 1.627 in. 

• The shaft diameter is 2.6798 in. 

• The reference clearance is 0.001 in. 

• The ratio of specific heats of the gas is 1.66 

• The gas constant is 1,790,000 in 2 /(s 2 -°R) 

• The absolute temperature is 528°R 

• The gas viscosity is 2.9 x 10" 9 lb-s/in 2 

• The shaft speed is 0 rpm 

• The reference pressure is 14.7 psig 

• The boundary pressures surrounding the seal are 50 psig 

e Cross -coupled stiffness and damping are to be computed at an excitation 
frequency of 0 rpm 

e There are 6 discrete inherently compensated source points in the sector 
of diameter - 0.020 in. The location of these orifices was determined from 
the design layout of the sector. The coefficient of discharge of each 
orifice is unity. 

e The buffer pressure is 200 psig 

e Flow is to be determined along four paths that make up the periphery of 
the seal . 

• The FILE option was exercised. A previous pressure distribution was read 
as the initial pressure distribution for this case. Convergence of the 
pressure is often difficult when solving source problems, whether they be 
inherently compensated sources or spot recesses. Convergence difficulties 
occur because sources present spikes in the pressure distribution and 
pressure gradients become very large. There are two methods for handling 
these problems which can be applied independently or jointly. The first is 
to use variable grid and fine grid spacing around the orifice holes. Each 
grid line around the hole should be at a distance of one to two orifice 
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SECTOR SEAL CROSS SECTION 
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Figure 5-20 T-Shaped Sectored Ring Seal 
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Axial Force Balance 
Fax “ Fp — ka 5 a = 0 


Figure 5-21 Pressure Distribution and Force Balance T-Sectored Seal 
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STiSIl 




diameters in both the axial and circumferential direction. The other 
mechanism is to start the problem at low eccentricity and use the pressure 
distribution as an initial guess to get to the next eccentricity. Continue 
the process until the desired eccentricity is attained. 

As indicated on Table 5-4, the eccentricity of the sector to support the load is 
0.55609 and the eccentricity angle is 90*. Figures 5-22 and 5-23 show the 
clearance and pressure distributions. Note on the plots the dense grid work 
surrounding the orifice locations. 
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TABLE 5-4 


Summary of Performance of T-Shaped Sectored Seal 
Sample Problem 4 


JOURNAL t LOAD POSITION 
ECCENTRICITY ■ 

ECCENTRICITY ANGLE » 

MINIMUM FILM « 

LOAD • 

LOAD ANGLE » 

POWER LOSS . ■ 

LEAKAGE AT I ■ 1 • 


LEAKAGE AT I ■ N 


STIFFNESS COEFFICIENTS 

PRINCIPAL X 

KXX ■ 

CROSS-COUPLED 

KXY ■ 

CROSS -COUPLEO 

ICXA > 

CROSS ‘COUPLED 

KXB ■ 

CROSS -COUPLED 

mrx ■ 

PRINCIPAL Y 

ICTY ■ 

CROSS -COUPLED 

KYA * 

CROSS -COUPLED 

KYB ■ 

CROSS -COUPLED 

KAX ■ 

CROSS-COUPLED 

ICAY « 

PRINCIPAL A 

KAA - 

CROSS -COUPLED 

KAB « 

CROSS -COUPLED 

ICBX - 

CROSS-COUPLED 

KBY * 

CROSS -COUPLED 

NBA « 

PRINCIPAL I 

KBB * 


.55609 
90.00 OEG 
.0004446 IN 
369.5 LB 

•90.00 OEG 

.0000 HP 

-.40231E-04 LB/S 

.40231E-04 LB/S 


.4316E+05 LB/IN 
.6956E-08 LB/IN 
•1299E-06 LB/RAO 
.2104E-02 LB/RAO 
•126.8 LB/IN 
.4982E+05 LB/ IN 
•295.4 LB/RAO 
•59.46 LB/RAO 
.3522E-06 1 N- LB/ 1 N 
-.3701E-02 IN-LB/IN 
.4361E+05 IN-LB/RAO 
-.3030E-06 IN-LB/RAD 
.1380E-02 IN-LB/IN 
-.9135E-10 IN-LB/IN 
.23271*08 IN-LB/RAO 
9965. IN-LB/RAD 


•DAMPING COEFFICIENTS 


PRINCIPAL X 

DXX ■ 

CROSS-COUPLED 

DXY ■ 

CROSS -COUPLED 

DXA ■ 

CROSS-COUPLED 

DXB ■ 

CROSS -COUPLED 

DYX « 

PRINCIPAL Y 

DYY « 

CROSS-COUPLED 

DYA » 

CROSS- COUPLED 

DYB * 

CROSS-COUPLED 

DAX * 

CROSS -COUPLED 

DAY » 

PRINCIPAL A 

DAA « 

CROSS-COUPLED 

DAB * 

CROSS-COUPLED 

D8X - 

CROSS -COUPLED 

DBY - 

CROSS -COUPLED 

DBA ■ 

PRINCIPAL B 

DBB * 

RIGHTING MOMENT 


ABOUT X-X 

NX ■ 

ABOUT Y-Y 

MY ■ 


10.82 LB-S/IN 
-.7990E-11 LB-S/IN 
-.2932E-12 LB-S/RAD 
•3324E-06 LB-S/RAD 
.1418E-01 LB-S/IN 
80.23 LB-S/IN 
.3588E-01 LB-S/RAD 
.6553E-02 LB-S/RAD 
.6908E-10 IN-LB-S/IN 
.8762E-06 IN-LB-S/IN 
4.516 IN-LB-S/RAD 
-.3370E-09 IN-LB-S/RAD 
-.2926E-06 IN-LB-S/IN 
•6752E-14 IN-LB-S/IN 
.4698E-12 IN-LB-S/RAD 
.9017 IN-LB-S/RAD 


.2753E-05 LB- IN 
-.6058E-14 LB-IN 


FLOW THRU SPECIFIED GRID LINE 

FROM 1 1 TO 27 1 FLOW* -.S908E-04 LB/S 


•FLOW THRU SPECIFIED GRID LINE 
FROM 1 34 TO 27 34 FLOW- .5903E-04 LB/S 


•FLOW THRU SPECIFIED GRID LINE 
FROM 1 1 TO 1 34 FLOW- -.4023E-04 LB/S 


-FLOW THRU SPECIFIED GRID LINE 
FROM 27 1 TO 27 34 FLOW- .4023E-04 LB/S 
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’CYL FILM THICKNESS DISTRIBUTION 

NASA SECTORED SEAL INHERENT COMPENSATION. DRAWING DIMEN 
DIAMETER = 2-680 IN SPEED = .00 RPM 

LENGTH = 1 -627 IN 



Figure 5-22 Clearance Distribution - T-Sectored Seal 
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Figure 5-2 3 Pressure Distribution - T-Sectored Seal 
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5,2^5 Sample Problem 5. Rav leieh-steo. Floatin g-Ring Seal 


This example represents another buffer fluid seal that was designed for use in 
the SSHE. The principal of operation of a hydrodynamic, lift-pad, floating-ring 
seal is illustrated on Figure 5-24. The seal consists of two rings that are 
mounted back-to-back. The buffer fluid enters between the rings and forces the 
rings up against the stationary housing. The buffer fluid leaks in the clearance 
annulus between the shaft and the seal and prevents ingress of exterior fluid on 
either side of the floating-ring assembly. The rings are held in equilibrium by 
a number of forces as shown on Figure 5-24. F c is a pressure force from the inlet 
buffer fluid that forces the rings up against the housings. This pressure force 
is partially balanced on the housing sides of the rings by undercutting and 
exposing the housing sides of the rings to buffer pressure. This balance force 
is identified as F B . F B represents a hydrodynamic force that is generated by 
rotation between the shaft and ring. The net hydrodynamic force is zero when the 
shaft and rings are in the concentric position. However, when the ring becomes 
eccentric with respect to the shaft, a hydrodynamic force is built up that 
opposes the eccentricity. There is also a normal force, F B , acting on the ring 
at the contact area between the ring and the housing. In addition to the 
equilibrium forces mentioned above, there is a friction force, F f , between the 
seal ring and housing. 

Figure 5-25 shows the hydrodynamic geometry that is incorporated into the bore 
of the seal rings. A portion of the length of the bore is segregated into 
sectors, and these sectors are separated from one another by axial grooves. A 
circumferential groove that goes completely around the bore is installed upstream 
of the final seal dam region. At the interior of the sectors, Rayleigh-step 
pockets are machined. The velocity direction of the shaft is such that it 
produces hydrodynamic pressures due to pumping of the fluid over the Rayleigh- 
step. The sealing occurs across the dam which is a narrow annulus of low 
clearance exposed to high pressure at its interior circumferential groove and to 
lower pressure at its outboard end. The shaded regions on Figure 5-25 indicate 
depressions from grooves and Rayleigh-step pockets. 

In this example one pad of the Rayleigh-step interface was examined from the high 
pressure interior end to the low pressure exterior end. The high pressure end 
is at the bottom end of the grid (I - 1). Geometric and operating parameters are 
as follows: 

e International units are to be employed 

e NPAD -1, because we are examining the performance of one pad only which 

is represented on the grid as a 90° arc from the center of one axial groove 
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Figure 5-24 Floating Ring, Rayleigh-Step Seal 
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Figure 5-25 Developed View of 50-mm Rayleigh-Step Pad 
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to the next. Thus the pad angle is 90°. The pad starts at 0°. Also, the 
boundary conditions at the circumferential ends of the single pad must be 
periodic, i.e. all pads will act identically, which will occur when the 
shaft is in the concentric position. Periodicity is invoked by applying 
the JOINED parameter. 

• The shaft diameter is .05 m 

• The total seal length is 0.0123 m 

• The seal clearance is 1.27 x 10' 3 m 

• The step height is 2.54 x 10 3 m 

• The gas viscosity is 2.19 x 10 3 N-s/m 2 

• The absolute temperature is 338. 6°K 

• The ratio of specific heats is 1.66 

• The gas constant Is 1154.8364 m 2 /(s 2 -°K) 

• The shaft speed is 70,000 rpm 

• The reference pressure is 101,352.93 Pa 

• The high pressure to be sealed is 1.37895 x 10 s Pa which would be 
at the bottom of the grid. The remaining boundaries are at 0 psig 

Table 5-5 shows performance of the sector examined. Since only one quarter of 
the seal is being examined, the leakages at the axial inlet <I-1) and outlet 
(I-M) would be multiplied by 4 to get total leakage. The leakage at the outlet 
is greater than the leakage at the inlet end because of the added flow 
contribution from the axial groove. The 3-D plots of the clearance and pressure 
distributions are shown on Figures 5-26 and 5-27 respectively. Note on Figure 
5-27 that the high ambient pressures in the grooves overshadows the increased 
pressure from the Rayleigh- step. 
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TABLE 5-5 


Performance Results 
Sample Problem 5 


•JOURNAL ft LOAD POSITION 
ECCENTRICITY * 

ECCENTRICITY ANGLE 
MINIMUM FILM * 

LOAD ■ 

LOAD ANGLE * 


.00000 
.00 OEG 
.0000127 M 
571.0 N 
-134.78 DEG 


POWER LOSS >21.79 U 


LEAKAGE AT I « 1 > -.21681E-04 KG/S 

LEAKAGE AT I > M - .40252E-03 KG/S 


•RIGHTING MOMENT 

ABOUT X-X MX « -.1596 M-N 

ABOUT Y-Y MY ■ .1551 M-N 
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FILM THICKNESS DISTRIBUTION 


GCYL 

RRYLEIGH - STEP SEAL PROBLEM 
DIAMETER = .050 M SPEED = 70000.00 RPM 


LENGTH = .012 M 

CLEARANCE = .000013 M 



Figure 5-26 Clearance Distribution - Rayleigh-Step Pad 
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GCYL PRESSURE DISTRIBUTION 

RAYLEIGH - STEP SEAL PROBLEM 


DIAMETER = .050 M SPEED = 70000-00 RPM 

LENGTH r .012 M 

CLEARANCE = .000013 M 



Figure 5-27 Pressure Distribution - Ray leigh-Step Pad 
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5.2.6 Sample Problem 6. Ravleigh-step S eal with Eccentricity 


This problem will be similar to problem 5 except the shaft is to be eccentric 
with respect to the seal ring. In this case periodic boundary conditions cannot 
be used, and to conserve grid space one hydrodynamic pad will be modeled and the 
number of pads will be four. To model separate pads however requires that the 
boundary conditions be known on all extremities of the pad. The seal dam region 
is not a separate pad problem but is a single 360° pad. Thus the problem resolves 
into two problems ; one that treats the separate Rayleigh pads and one that treats 
the seal dam. For this particular example, only the Rayleigh- step hydrodynamic 
region is considered. The following geometric and operating parameters have been 
applied: 

• International units are invoked 

e OPTION -1, the shaft position relative to the seal ring is specified 

• Stiffness is to be calculated in four degrees of freedom at an 
excitation frequency of 70,000 rpm. 

e The number of pads is 4 and each pad has an extent of 90° 
e The shaft diameter is 0.05 m 

• The shaft length is 0.0123 m 

e The reference clearance is 1.27 x 10 5 m 

• The gas viscosity is 2.19 x lO'^-s/m 2 
e The absolute temperature is 338. 6°K 

e The ratio of specific heats is 1.66 
e The gas constant is 1154.84 m 2 /(s 2 -°K) 

• The shaft eccentricity ratio is 0.5 
e The eccentricity angle is 270° 

e The shaft speed is 70,000 rpm 

• The reference pressure is 1.01353 x 10 3 Pa 
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• The pad boundary pressures are 1.37895 x 10 6 Pa 


As shown on Table 5-6, at the specified position, the load capacity of the seal 
is 52.96 N and the load angle is 71.63° from the x-axis. The minimum film 
thickness is 6.4 x 10' 6 m. The clearance and pressure distributions are shown 
on Figures 5-28 and 5-29 respectively. 
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TABLE 5-6 


Summary of Results 
Sample Problem 6 


-JOURNAL C LOAD POSITION 
ECCENTRICITY ' 

ECCENTRICITY ANGIE 
MINIMUM FILM 
LOAD 

LOAD ANGLE 
POWER LOSS 
LEAKAGE AT I > 1 
LEAKAGE AT I ■ H 
•STIFFNESS COEFFICIENTS 


PRINCIPAL X KXX 

CROSS -COUPLED KXY 

CROSS-COUPLED KXA 

CROSS-COUPLED KXB 

CROSS-COUPLED KYX 

PRINCIPAL T KYT 

CROSS-COUPLED KYA 

CROSS-COUPLED KYB 

CROSS-COUPLED KAX 

CROSS-COUPLED KAY 

PRINCIPAL A KAA 

CROSS-COUPLED KAB 

CROSS-COUPLED KBX 

CROSS-COUPLED KBY 

CROSS-COUPLED KBA 

PRINCIPAL B KBB 

•DAMPING COEFFICIENTS 
PRINCIPAL X DXX 

CROSS-COUPLED DXY 

CROSS-COUPLED DXA 

CROSS-COUPLED DXB 

CROSS-COUPLED DYX 

PRINCIPAL Y DYY 

CROSS-COUPLED DYA 

CROSS-COUPLED DVB 

CROSS-COUPLED DAX 

CROSS-COUPLED DAY 

PRINCIPAL A DAA 

CROSS-COUPLED DAB 

CROSS-COUPLED DBX 

CROSS-COUPLED DBY 

CROSS-COUPLED DBA 

PRINCIPAL B DBS 

•RIGHTING MOMENT 
ABOUT X-X MX 

ABOUT Y-Y MY 


.50000 
-90.00 0E6 
.0000064 N 
52.96 N 

71.63 DEG 

95.08 W 

-.88483E-04 KG/S 

.88483E-04 KG/S 


.1016E+08 N/H 
•3371E+07 N/N 
14.33 N/RAD 
-.3404 M/RAD 
-.6161E+06 N/N 
.1169E+08 N/N 
37.51 N/RAD 
20.59 N/RAD 
-.8526E-10 N-M/M 
-.S781E-10 N-M/M 
27.53 N -M/RAD 
1.657 N-N/RAD 
.1284E-09 N-M/M 
-.8787E-11 N-M/M 


*10.32 

N-M/RAD 

14.46 

N-M/RAD 

786.4 

N-S/M 

*216.6 

N-S/M 

-.4094E-03 N-S/M 


.8221E-04 N- S/RAD 

192.7 N-S/M 

901.7 N-S/M 
-.101 IE-02 N-S/RAD 
-.6955E-03 N-S/RAD 

.4642E-15 N-N-S/N 
.5113E-15 N-N-S/N 
.3S72E-02 N-M-S/RAD 
-.1573E-03 N-M-S/RAD 
-.9665E-15 N-M-S/M 
.2712E-15 M-M-S/M 
-.5069E-04 N-M-S/RAD 
.2612E-02 N-M-S/RAD 


.7687E-15 N-M 
.4289E-15 N-M 
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JCYL FILM THICKNESS DISTRIBUTION 

RAYLEIGH - STEP SEAL PROBLEM 

DIAMETER = .050 M SPEED = 70000-00 RPM 

LENGTH = .012 M 

CLEARANCE = .000013 M 



Figure 5-28 Clearance Distribution - Rayleigh-Step Seal with Eccentricity 
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Figure 5-29 Pressure Distribution - Ray leigh-Step Seal with Eccentricity 
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5.3 Verification of GCYL 


Several mechanisms were used to conduct verification of the code. For the most 
part the results of the code were compared against information in the public 
domain literature, and in some instances, comparisons were made against the 
results of other codes and against manual computations. 

The first case is the pressure distribution of an infinitely long slider. The 
results are compared at several values of A as shown on Figure 5-30. 

Further comparisons were made for a plain cylindrical seal with an L/D ratio of 
1 with information from reference 3. Computations were made at two different 
eccentricity ratios, € - 0.6 and 0.8. Non dimensional load capacity and attitude 
angles are shown on Figures 5-31 and 5-32 respectively. Excellent correlation is 
demonstrated. 

A significant feature of the GCYL code is the computation of frequency dependent 
stiffness and damping coefficients. The method was first implemented in the 
compressible Spiral-Groove computer code SP1RAIX?. These stiffness and damping 
coefficients are important because they are used to represent the fluid film 
characteristics in dynamic analysis . Their computation embodies many features of 
the code including steady state performance. Table 5-7 shows comparisons for 
three codes for an excitation frequency of zero, and for a 360° cylindrical seal 
in the concentric position. The first column represents the code GCYL as 
previously modified with only the capability to compute zero excitation 
frequencies. The second column represents the latest version of the code with the 
frequency dependent stiffness and damping routines. The third column are the 
results produced by the Spiral Groove code with zero groove depth, so that the 
geometry of the three cases are equivalent. 

Table 5-7 

Comparison of Spring and Damping Coefficients 


Coefficient 

Unit 

GCYL 

(Previous) 

GCYL 

(New) 

SPXRALG 


lbs/in 

5,715 

5,719 

5,715 


lbs/in 

7,301 

7,301 

7,363 


lbs/ in 

-7,107 

-7,092 

-7,140 

*yy 

lbs/in 

12,550 

12,590 

12,752 

K.. 

in- lbs/rad 

384.2 

384.2 

395 
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0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 


x/L 


Figure 5-30 Rayleigh-Step, Program Verification 

*"Theory of Hydrodynamic Lubrication", 0. Pinkus, B. Stemlicht, 
McGraw-Hill, New York, 1961 
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Plain Cylindrical Seal, L/D = 



(aid)/M 
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Plain Cylindrical Seal, L/D = 



<d‘a|6uv apnjjnv 
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Design of Gas Bearings", Published by Mechanical Technology Inc 






Comparison of 

Table 5-7 

Spring and Damping Coefficients 


Coefficient Unit 


GCYL 

(Previous) 

GCYL 

(New) 

SPIRALG 


in- lbs/rad 


181.6 

181.6 

194 

Kt. 

in- lbs/rad 


-277.2 

- 277.2 

-300 

Kbb 

in- lbs/rad 


122.6 

122.6 

124 

Damping Constants 


lb-s/in 


1.402 

1.387 

1.403 


lb-s/in 


-1.888 

-1.870 

-1.855 

Dyx 

lb-s/in 


2.992 

2.999 

3.019 


lb-s/in 


1.957 

1.897 

1.901 

D.. 

in- lb -s/rad 


0.1088 

0.1079 

0.1148 

D* 

in- lb- s/rad 


-.0447 

-.0450 

-.046 

Db. 

in - lb -s /rad 


0.0377 

0.0377 

0.0385 

Dbb 

in- lb- s/rad 


0.0662 

0.0659 

0.0695 

L - 1 in 
V 0.5, 

. , D - 1 in. , C - .001 
excitation frequency — 

in. , /i - 3x10"’ 
0 

lb-s/in 2 , N - 

48,000 rpm 


Kjj is the stiffness in the i direction due to a j displacement 
is the damping in the i direction due to a j velocity 
x and y are translations and a and b are rotations. 

Table 5-8 shows a comparison for a synchronous excitation (excitation frequency 
that is equal to the shaft speed). For both situations, the correlation is very 
good. 
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Table 5-8 

Stiffness and Damping Comparison at Synchronous 


Coefficient Unit GCYL 


Frequency 


Unit 

GCYL 

SPIRALG 

lbs/in 

7,498 

7,467 

lbs/in 

1,122 

1,203 

lbs/in 

-1,122 

-1,203 

lbs/in 

-7,498 

-7,467 

in- lbs/rad 

100.1 

99 

in- lbs/rad 

127.7 

134 

in- lbs/rad 

-127.7 

-134 

in- lbs/rad 

100.1 

99 



Table 5-9 shows the variations in stiffness and damping for a 360° cylindrical 
seal with an excitation frequency equal to operating speed as compared to an 
excitation frequency of zero. 
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Table 5-9 

Stiffness and Damping Coefficients at 
Two Excitations 


Operating Speed 

rpm 

48,000 

48,000 

Excitation 

rpm 

48,000 

0 

Frequency 




K** 

lbs/in 

9,648 

5,885 

K*y 

lbs/in 

1,942 

7,267 

K« 

lb/rad 

0.7303 

-0.2391 


lb/rad 

1.291 

1.298 

Kyx 

lb/in 

1,040 

-7,116 

Kyy 

lb/in 

17,670 

13,050 

Ky. 

lb/rad 

1.451 

3.258 

Kyb 

lb/rad 

0.4477 

1.877 

K.. 

in- lb/rad 

639.4 

420.6 


in- lb/rad 

71.58 

192 

Kb. 

in- lb/rad 

-193.8 

-293.8 

Kbb 

in- lb/rad 

221.1 

133.7 

Damping Coefficients 

Dxx 

lb -s/in 

1.658 

1.406 


lb -s/in 

-0.7059 

-1.859 

D y* 

lb -s/in 

0.9180 

3.012 

Dyy 

lb -s/in 

1.521 

1.897 

D.. 

in- lb- s/rad 

0.090 

0.113 


in- lb -s/rad 

-0.0311 

-0.0473 

Db. 

in- lb -s/rad 

0.0256 

0.0393 

°bb 

in- lb -s/rad 

0.0666 

0.069 

360° Cylindrical Seal, L - 1 in, D - 
H - 3 x 10' 9 lb- s/in 2 , 0 gage pressure 

1 in, C - .001 in, 
at both ends 
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Substantial differences are noted, which demonstrates the significance of 
applying proper frequencies when computing stiffness and damping. 


An internal check of the code can be made by analyzing a recessed hydrostatic 
bearing. With the flow path option, the net flow around the periphery of a 
hydrostatic pad can be determined and compared against the inflow to the recess. 
For flow continuity, the sum of the peripheral flows should equal the inlet flow. 
The following geometry and operating parameters were considered. 

e A single pad with grid dimensions of 15 x 37 (M x N) . 
e The pad diameter is 2 inches 
e The pad length is 2 inches 

• The pad clearance is 0.001 in. 

• The pad angle is 180° and the starting angle is at 180° 

• There is one recess located in the pad, and the grid corner points are 
as follows: Left bottom corner, M - 3, N - 22 

Right Top corner, M - 13, N - 27 

• The specific heat of the gas is 1.4 

• The gas constant is 250,000 in 2 /(s 2 -°R) 

• The absolute temperature is 530° R 

e The absolute viscosity is 3 x 10" 9 lb-s/in 2 

e The inlet orifice diameter to the recess is 0.020 in and the coefficient 
of discharge is 1.0. The orifice is located in the grid at M - 8, N - 24. 

• The supply pressure to the orifice is 150 psig. The pressure surrounding 
the pad is at 0 psig. The reference ambient pressure is 14.7 psia. 

• Several eccentricities and speeds were examined and are defined in the 
subsequent discussions. 

The output from the code supplies the total flow from the peripheral flow path 
and the pressure in the recess. A manual computation can then be made for 
calculating the inlet flow through the orifice using the following equation: 


f 0 = 386.4 AoC&p, 


l l p *J i ) 


(5-88) 


where , 


I 


p. — 

G c 6(Y-1) 


(5-89) 
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f 0 - inlet flow, lb/s 
A„ - orifice area, in 2 
C D - discharge coefficient 
p, - supply pressure, psia 
p r - recess pressure, psia 
y — ratio of specific heats 
G c “ gas constant, in 2 /(s 2 -°R) 
0 - absolute temperature, °R 


Table 5-10 provides the 

results of several 

cases 




Table 

5-10 Recessed Pad Flow Comparisons 


€ 

N 

Q P 

Pr 

Qo 

A 


rpm 

lbs/s 

psig 

lbs/s 

X 

0.0 

0.0 

0.001188 

43.2* 

0.001189 

0.08 

0.4 

0.0 

0.001174 

84.4 

0.001176 

0.17 

0.0 

70,000 

0.001188 

36.1* 

0.001189 

0.08 


* Choked Flow 

€ — Eccentricity ratio 

N — Shaft speed 

Q p - Peripheral flow 

p r — Recess pressure 

Q,, .Orifice flow 

A - percent variation 

Note that the peripheral and orifice flows differ by less than 0.2Z. 

When using the source points or spot recess options of the code, it is important 
to surround the source point with a fine grid to obtain an accurate result and 
a computation in which pressures will converge. Studies were made of varying grid 
sizes for a source problem. The variable grid option was applied and varied. A 
single pad with a central row of orifices were analyzed (see sample problem 
Number 4). The following information is pertinent: 

• Number of pads - 1 

• Pad angle - 120° 

• Start angle - 30° 

• Number of grid points in circumferential direction - 37 

• Number of grid points in axial direction - 15 
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• Diameter — 2.6798 in. 

• Length - 1.627 in. 

• Specific Heat Ratio - 1.66 

• Gas constant — 1,790,000 in 2 /( s2 *°^) 

• Absolute Temperature - 528°R 

• Viscosity - 2.9 x 10 '* lb -s/in 2 

• Shaft Speed - 0 rpm 

• Reference pressure — 14.7 psia 

• Boundary pressures - 60 psig 

• Supply pressure to inherently compensated orifices 

• Preload - 20X located at the center of the pad 

• Stiffness is to be determined 

• Six source points are located along a circumferential line in the axial 
center of the pad at circumferential grid locations 5, 10, 15, 20, 25, 30. 

The hole diameter is 0.015, and the coefficient of discharge is 1.0. 

Tables 5-11 and 5-12 indicate the effect of grid width around the source point 
in both the axial and circumferential directions. Table 5-11 indicates the 
source pressures as the grid width is changed. They are relatively unaffected 
until the grid width is 8x the orifice hole size. A similar conclusion can be 
drawn for the other performance parameters of load, flow, stiffness and damping 
as indicated by Table 5-12. The recommended grid width from the source point to 
a neighboring grid line is twice the orifice diameter. 
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Table 5*11 

Comparative Studies - Discrete Orifices Vs. Grid Size 
Orifice Size - 0.015 in. 

A — Grid width around orifice in both circumferential and axial directions 

Comparison of Source Pressures 


A 

in 

Pi 

psig 

P 2 

psig 

P 3 

psig 

P* 

psig 

P 5 

psig 

P* 

psig 

0.015 

129 

144 

151 

151 

145 

129 

0.030 

128 

144 

151 

151 

144 

128 

0.060 

126 

143 

149 

149 

143 

126 

0.120 

122 

140 

147 

136 

140 

122 


Table 5-12, Comparison of Performance 


A 

in. 

U 

lbs 

Qi 

lbs/s 

Qm 

lbs/s 

*xx 

Lbs/in 
x 10' 6 

lbs/in 
x 10' 6 

Dxx 

(lbs- 

s)/in 

Dyy 

(lbs- 

s)/in 

0.015 

314.9 

0.11338 

0.11338 

0.0419 

0.1320 

2.945 

16.43 

0.030 

315.4 

0.11393 

0.11393 

0.424 

0.1338 

2.933 

16.33 

0.060 

316.4 

0.11496 

0.11496 

0.0432 

0.1371 

2.911 

15.83 

0.120 

319.3 

0.1186 

0.1186 

0.0446 

0.1480 

2.866 

11.75 


A - grid width 
W - load capacity 
Q x — flow out of grid line M-l 
flow out of grid line M-M 

and Kyy - Stiffness in x and y directions, respectively 
and Dyy - Damping in x and y directions, respectively 
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5.5 NOMENCLATURE 


C(j = inherently compensated orifice coefficient of discharge 
C Q - reference clearance (concentric clearance) 
d 0 = orifice diameter 

e - shaft displacement from concentric position 

Ff * viscous friction force 

FF * dimensionless viscous friction force = Ff/(p 0 C 0 R) 

Gq = universal gas constant 

h s local film thickness 

H * dimensionless film thickness = h/C 0 
1 * bearing length 

L * dimensionless length = 1/R 

N * number of orifices in a row 


P 

Po 

P 

p CR 

PR 

p S 

q 

R 

r 

S c 

t 

t 0 


pressure 

reference pressure 

dimensionless pressure = p/p Q 

critical pressure ratio 

orifice downstream pressure 

supply pressure upstream of orifice 

mass flow 

journal radius 

orifice hole radius 

source correction factor 


time 

reference time * 



T s dimensionless time * t/t 0 

T a = absolute temperature 

Tf = viscous friction torque 

TF = dimensionless viscous friction torque * Tf/(p 0 C 0 R^) 
U = journal surface velocity 

z » axial direction coordinate 

Z = dimensionless axial coordinate = z/R 
a s misalignment angle about x-x axis 

B = misalignment angle about y-y axis 

Y = ratio of specific heats 
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e 


e 

e 

p 

A 


= eccentricity ratio = e/C 0 

= angular direction (direction of sliding) 

= angular extent of pad 2 

6ptoR 

= compressibility parameter = 2 

Po^o 


M = absolute viscosity 

€a> = rotating speed 
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6.0 Knowledge Base System Development 

One of the significant aspects of the overall program is the generation of a Knowledge Based 
System (KBS) having the following objectives: 


1) To integrate the scientific and industrial codes into a package that will 
provide access to important technical data and information to facilitate 
generation of optimum seal configurations. 

2) To provide a user friendly graphical user interface with context sensitive 
help. 


3) To provide Expert systems to help select the type of seal best suited for the 

intended application, analyze user input and output of analysis codes to 
guide the seal design optimization process. 


This report describes the architecture of the KBS and the development of the user interface 
elements during the first year of the program. 


6.1 KBS COMPONENTS 

A schematic of the KBS is shown in Figure 6-1. Functions of the various components are 
described below. 


6.1.1 Executive Program 

The executive shell integrates all the components of the KBS and provides the user with a 
single point of access for all the resources in the KBS. Features of file executive are: 


■ Access to scientific and industrial codes. 

■ Access to the expert systems for seal type selection and seal design 
guidance. 

■ Utility functions including browsing and printing output files created by 
the analysis programs, plotting routines to display the results in a graphical 
form, and procedures for the users to add their own programs to the KBS. 

■ Network communications with the Cray X-MP computer used for running 
the scientific codes. The communication procedures will be made as 
transparent as possible. 



NASA Workstation 
Functional Block Diagram 



Scientific Scientific 

Code Code 











■ Database services to access the databases vised to store input and output 
data sets for the analytical codes. The access will be controlled using 
passwords to prevent unauthorized access. 


6.1.2 Scientific Codes 

Scientific seal analysis codes will provide steady-state and transient analysis capability based 
on full three-dimensional Navier-Stokes equations. Other characteristics of these codes 
include the following: 

■ Cylindrical, polar, and non-orthogonal body-fitted coordinates 

■ Stationary and rotating coordinate systems 

■ Advanced turbulence models suitable for high-shear rotating flows 

■ Incompressible and compressible flows 

■ Cavitation or liquid film rupture 

■ Energy conservation equation with viscous heating and phase changes 

■ Provisions for additional field equations such as electromagnetic and 
electrostatic forces 

The seal a nal ysis will encompass a number of generic seal categories including Cylindncal, 
Labyrinth, Damper, Honeycomb, Face, Noncontinuous, Wave, Grooved, Tip, Contact, and 
Brush seals. The models for these seals will be very extensive and detailed, requiring a Cray 
X-MP class computer for execution in a reasonable time. 

The KBS will be used to prepare the input data and the input files sent to the Cray using a 
network. After execution, the output files will be downloaded from the Cray computer and 
post-processed on the KBS. The users of the scientific codes are expected to be high-level 
research personnel familiar with the fundamental theories used in five codes, the 
mathematical underpinnings and the basic structure of the code, and the types of seals being 
analyzed. The assistance provided to these users will focus on the mechanics of defining the 
analytical model during input. An extensible database of typical models will be provided as a 
starting point for user input. The final user input will be checked using expert systems to 
ensure that all the necessary input has been supplied and that there are no obvious errors in 
grid specification, boundary conditions, material characteristics, etc. This function is essential 
given the expense of running these codes on a Cray computer. Expert systems and 
conventional data reduction software will be provided to assist the user in interpreting the 
outputdata. 
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6.1.3 Industrial Codes 

The industrial codes are simpler two- and three-dimensional codes for several different types 
of seals. Some of the codes included in the package are: 

■ Bushing and Ring Seals 

❖ Uniform 

❖ Axial Step and Taper 

❖ Hydrodynamic Step and Taper 

❖ Self Energized Hydrostatic 
«5* Segmented 

■ Face Seals 

❖ Contact Face Seals 

❖ Radial Step and Taper 

❖ Hydrodynamic Step and Taper 

❖ Hydrostatic 

❖ Spiral Groove 

❖ Multi-pad 

■ Labyrinth Seals 

❖ Straight 

❖ Stepped 

❖ Abradable 

❖ Angled 

■ Tip Seals 

■ Damping Seals 

■ Brush Seals 

■ Electro-fluid Seals 


Smart Seals 



The anticipated users of the industrial codes include seal design and application engineers, 
seal users such as rotating machinery designers, and analysts performing seal design audits 
and failure analysis. These users may not be familiar with all of the seal analysis capability 
incorporated into these codes or with all of the design and analysis options available for a 
specific application. The assistance provided to these users will include extensive on-line, 
context-sensitive help for each code, and error trapping to ensure that input values are within 
admissible limits. A graphical user interface using windows and drop-down menus will be 
provided for each code to ensure a uniform look and feel. The names of menu items and seal 
variables will be standardized to reduce the learning curve. Expert systems will be provided 
as needed to guide the user in setting up an optimum analytical model and to interpret the 

output data. 

Seal design optimization is an iterative process involving seal interfacial analysis, 
rotordynamic analysis, and thermo-elastic analysis. Expert systems will be provided to guide 
the user through this process by helping the user to select the analyses to be performed and 
to decide when to terminate the iterative process. 


6.1.4 Databases 

Databases will be included to store input and output data sets for example problems and for 
problems used for analytical code validation. The databases will also enable users to develop 
a library of analytical models tailored for their individual needs and to maintain a history of 
analyses performed using the codes. The analytical codes will access their databases directly 
to store and retrieve data. Users may also access the databases using database services 
provided through the executive. 

6.2 HARDWARE AND SOFTWARE SYSTEM SELECTION 
6.2.1 Initial System Selection 

The hardware and software system selection were driven by the need for a graphical user 
interface and the computational requirements for the analysis codes. The choice was 
complicated by having to anticipate the hardware and software availability six years down 
the road, by a wide variation in the computing power available to the anticipated users, and 
by the wide range of computational requirements for the individual analysis codes. 

After a review of KBS requirements, an Intel 80386 and 80486 based IBM PC compatible 
hardware platform running OS/2 with the Presentation Manager interface was selected. 
FORTRAN 77 was selected for implementing the analytical codes and C as the primary 
language for user interface development. Developing graphical user interfaces is a time 
consuming process. The Toolbook authoring system was selected to explore alternatives to 
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developing the user interfaces in C. NEXPERT Object was selected as the expert system shell. 
NEXPERT is available for several operating environments and provides portable knowledge 
bases. Other options considered included UNIX on RISC workstations and Windows 3.0. The 
reasons for selecting OS/2 were as follows: 


■ The ease-of-use features planned for the KBS provide the most benefit for 
projected users of the industrial codes who may not be intimately familiar 
with the content of the codes and may not be comfortable using computers. 
These users are in organizations that typically use Intel based, IBM PC 
compatible machines. Therefore, Intel 80386 or 80486 based machines were 
selected to enable users to use existing hardware. 

H Windows environment was not acceptable because of its DOS based 
limitations on memory, networking support, etc Most of the analytical 
codes require more resources than provided by DOS. While it was possible 
to provide a collection of utilities to overcome these shortcomings, the cost 
of maintaining and developing a large software package like the KBS with 
such makeshift arrangements would be prohibitive in the long run. 

■ UNIX is a complex multi-user, multitasking operating system. It requires 
considerable expertise to install and manage UNIX systems. The expected 
users of the industrial codes are used to simple DOS systems and usually 
do not have the expertise or the support staff to manage UNIX systems. The 
multi-user capabilities provided by UNIX do not add enough value to 
offset the added complexity and cost of development. 

■ OS /2 combines the best features of UNIX and Windows environments. 

■ OS/2 is a single-user system. This reduces the complexity of die operating 
environment and the time required to learn it. However, the system 
security features such as log in control provided by UNIX and needed for a 
multi-user system are not available in OS/2. 

■ Like UNIX, OS/2 is a robust multitasking system designed for use over 
networks. In addition, it has features such as Multi-threading, Dynamic 
Data Exchange (DDE), Dynamic Link Libraries PLLs), and Installable File 
System (IFS) which are not yet available in UNIX. These features implement 
operating system concepts developed in recent years and provide a flexible, 
extensible operating environment. Multi-threading is essential for support 
of multi-processor hardware platforms. DDE allows for data links between 
programs that are more flexible than traditional inter-process 
communication method. DLLs allow development of code modules shared 
by several programs and easy upgrade of software already in use. 

■ The cost of third party software for OS/2 is comparable to DOS and 
Windows environments, and is much cheaper than similar UNIX software. 
The ability of OS/2 to run current DOS software in the DOS compatibility 
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mode and the use of the DOS file system makes for an easy transition for 
expected users of the industrial codes. 

■ OS/2 Extended Edition is the most comprehensive software development 
environment available today. The Application Programming Interfaces 
(API) for the base operating system, the database manager, and the 
communications manager provide a rich, robust, and well integrated 
environment to develop applications with seamless access to databases and 
networking capabilities. On any other platform these capabilities would 
have to be duplicated with tools from several different vendors which do 
not always work together cleanly. That is why OS/2 is becoming the 
platform of choice for mission critical applications and for downsizing 
mainframe or minicomputer applications using a client-server architecture. 

■ OS/2 Version 2.0, due before the end of 1991, will be a 32-bit operating 
system with an API that is portable to non-Intel hardware platforms. 


Based on comments by the peer review panel after the first workshop at NASA Lewis 
research Center, development of the KBS was interrupted to re-evaluate the harware and 
software system selection. A segment of the attendees at the annual workshop and members 
of the Peer Review Panel appointed by NASA favored the UNIX operating system. NASA 
did not want to deny UNIX users access to the KBS by limiting the development to OS/2. 


6.2.2 Hardware and Software System Re-evaluation 

Issues of software portability, cost of a delivery system, availability of adequate software 
development and maintenance tools, cost of software development, and expected evolution 
in operating system environments are under evaluation. 


6.2.2. 2 Hardware Platforms 

The hardware choices being evaluated are Intel 80386 or 80486 based IBM PC compatibles 
and RISC workstations. The advantages of an Intel-based platform are wide availability and 
low cost. Most of the anticipated industrial users are already using Intel-based machines. The 
major advantage of RISC machines is better floating-point performance. There are, however, 
several competing RISC platforms. This makes it necessary to commit to one hardware 
vendor or bear the additional costs of supporting all the different platforms used by the 
anticipated users of the analysis codes. It is anticipated that the performance of high-end PCs 
and low-end workstations will converge over the next few years. Therefore, the hardware 
platform should not be a decisive factor in the final system selection process. 
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6.2.22 Operating Systems 

The two operating systems being evaluated are a flavor of UNIX with OSF/MOTIF user 
interface, and OS/2 with the Presentation Manager interface. Both these systems offer a 
multi-tasking environment which is a must for a system that has analytical codes that take a 
long time to execute. The user must have access to other functions on the machine while 
analysis proceeds in the background. 


The advantages of UNIX are portability between Intel-based and RISC platforms running the 
same version of UNIX. Major disadvantages of UNIX are the number of different versions 
available from different vendors and the high cost of third-party software. UNIX vendors 
seem to be converging on two versions: System V Release 4 from Unix International with the 
Open Look Interface and OSF Unix with the Motif interface. Porting software between 
versions is not a trivial task. Additional disadvantages of UNIX on Intel-based machines are 
slow speed, large memory and disc size requirements. 


The advantages of OS/2 are a modem, integrated design that provides all the facilities 
needed for developing applications with graphical user interfaces, ability to run thousands of 
existing PC-DOS applications, lower cost third party software, and lower memory and disc 
size requirements. Currently, the major disadvantage is lack of portability to non-Intel 
hardware. However, a portable version of OS/2 that uses the same programming model as 
OS/2 version 2.0 is expected to be available in the 1992-93 time frame. 

Apple Computer and IBM recently signed a letter of intent to jointly develop a new 
object-oriented, portable operating system. A new object-oriented operating system that will 
run on Intel x86. Motorola 680x0, and IBM RS/6000 based computers and which will run 
existing ADC, Macintosh, and OS/2 applications. The new platform will be developed by a 
new company jointly owned by Apple and IBM. An industry standards group will be formed 
to set standards for the emerging software architecture and act as a clearinghouse for 
information about the project The operating system will be available in two to three years 
and will be made available to other hardware platform vendors. The operating system will be 
based on reusable software components easily portable to various systems. The object 
oriented nature of the system will allow vendors to differentiate their version of the 
enviroments by adding features such as portions of Macintosh and OS/2 operating systems 
or other applications. The level of customization will be greater than that provided in the 
current NextStep environment on NEXT machines, which is the only object-based operating 
environment currently available. Object oriented technologies developed by the joint 
company will be incorporated into OS/2 and the Macintosh operating systems as it becomes 
available to facilitate their integration into the jointly developed operating system. Apple and 
IBM will provide application programming interfaces (API) for the older operating systems 
like OS /2 and Macintosh to allow them to run new object-oriented applications. For IBM, 


6-8 


OS/2 will evolve to be a migration path to the new operating system in the late 1990’s. An 
encapsulation technology will be provided to allow older applications to run on the new 
hardware/software platforms. This may include providing full binary compatibility so that 
the older applications do not have to be recompiled. 


6.2.23 Portability Issues 

Software portability was considered at three different levels: 


1. The analysis codes are being developed using ANSI standard FORTRAN 77 
for portability. The codes are currently being designed to read input from a 
file and write output to a file. These codes need only to be recompiled for 
use on any system. 

2. Input and output post-processing programs with an easy-to-use graphical 
user interface. These programs will create file input files needed by an 
analysis code and read the output files for post-processing, viewing results, 
etc. This portion of the code is usually operating system dependent. 
Software tools that provide portable code for porting graphical user 
interfaces between UNIX and OS/2 are now becoming available. These are 
being considered. 

3. Advanced, interactive two- and three-dimensional graphics capability. This 
capability is envisioned for future versions of file analysis codes as 
low-cost, portable software tools become available for developing 
interactive graphics applications. The current graphical user interfaces are 
being designed to enable future incorporation of these capabilities. 


6.2.23 Software Development Tools 

MTI experience with developing the user interface for the industrial codes showed that 
conventional approaches to developing software are not cost effective when developing 
software with graphical user interfaces. A number of tools are becoming available that 
significantly reduce the cost of developing a graphical user interface and, in some cases, 
provide means to port the interface between different operating environments such as OS/2 
and UNIX with OSF/Motif interface. These packages usualy allow porting of user interface 
elements such as windows, menus, dialog boxes, scroll bars, and buttons. The code for other 
graphics elements drawn using graphics primitives such as lines, areas, curves, and special 
features such as detectable, dynamic segments for interactive graphics is not ported. These 
elements can be implemented using function libraries that provide support for standards like 
PHIGS and GKS. However, these libraries usually require a run-time license for each 
workstation that runs code using the library. Some of the available tools for developing 
graphical user interfaces are described below. 
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Open Interface from Neuron Data. Open Interface is a software development environment 
that allows the development of portable code for a graphical user interface for 
DOS/ Windows, OS/2 Presentation Manager, Unix with Motif and Open Look, Macintosh, 
and VAX VMS DEC Windows environments. The interface code is developed using a 
graphical screen layout tool which generates codes using generic C function calls designed 
by Neuron Data. This code is then compiled and linked with libraries for die environment to 
which the code is being ported. Cost of the development system is $9,000 for OS/2 and 
$12,000 for Unix. Run-time systems are $350 and $500 for OS/2 and Unix, respectively. The 
Open Interface environment was used by Neuron Data to develop the NEXPERT Object 
expert system shell that is availble for all of the above environments. 


CaseWorks. CaseWorks is a software development environment for developing a graphical 
user interface for Windows 3.0, OS/2 Presentation Manager and Unix Motif environments. 
The development system generates C code for an interface defined using CaseWorks tools. 
The Unix Motif product has been developed but is not shipping at this time because of a lack 
of demand for it. The company claims it is ready to ship it if there is sufficient demand for it. 
The cost of the OS/2 version is $2,000. There are no run-time fees. 


Extensible Virtual Toolkit (XVT). XVT is a set of libraries that support the Macintosh, 
Windows, OS/2 Presentation Manager, Open Look, and Motif environments. The user 
interface is developed using XVT function calls, and the resulting code is copiled and linked 
using an XVT library for the environment in which the application is going to run. The cost 
of XVT is $800 for DOS, Macintosh and OS/2 environments and $3,500 for UNIX 
environments. 


Information Engineering Facility (IEF) from Texas Instruments. IEF is a DOS and OS/2 
based CASE tool. The current version supports software development for OS/2 Presenta on 
Manager and Windows graphical user interfaces. A new version due out at the end of the 
year will let the developers distribute applications developed with IEF on HP 9000 and IBM 
RS /6000 running UNIX with the MOTIF GUI. This facility is geared towards developing 
business software that uses a mainframe database as a central source of data, with software 
running on DOS and OS/2 workstations. Prices vary depending on the configuration but are 

in the $6,000 to $10,000 range. 


Gpf - GUI Programming Facility from Microformatic. Gpf is a software development 
environment for developing a graphical user interface for the OS/2 Presentation Manager. 
The development system generates C code for an interface defined using Gpf tools. Gpf is a 
sophisticated tool that includes tools for reading data from the OS/2 Database Manager 
databases. Gpf generates the code with the SQL commands and the C function calls required 
to open and read data from the databases. The cost of Gpf is $3,500 with no run-time fees. 


6-10 

cV 



Gpf is not available for any other platform at this time. It would be the tool to use if we stay 
with OS/2 as the operating system. 


Expert System Shells. The CUPS library of C routines available from COSMIC is currently 
limited to a forward-chaining reasoning capability. An object-oriented version is being 
developed. However, CUPS is best suited for embedded expert systems and may not be 
convinient for this project if expert system capabilities are to be integrated with advanced 3-D 
interactive graphics in future versions of the codes. We need an object oriented expert system 
shell that is potable across several operating environments. NEXPERT Object seems to be the 
best compromise based on capabilities, portability, and cost Hurd party tools are available 
for developing interactive graphics applications using NEXPERT. The selection of an expert 
system shell can be put off until early 1992. 

Database Management Systems. The cost of OS/2 Extended Services (this includes the 
Database Manager, the Communications Manager and the OS/2 LAN Requester) is about 
$600. IBM has announced that the OS/2 Database Manager is being ported to AIX and will 
be available in the second or third quarter of 1992. The price for the ADC version has not been 
announced. Third party relational database management systems such as Oracle and 
Informix are available for both OS/2 and UNIX. Cost of a single user version of Oracle is 
$2,000 for UNIX and $1,500 for OS/2. A C language interface is included in the price. The 
run-time version of Oracle for OS/2 costs $200; The price for a UNIX run-time version was 
not available. The cost for the required took for a development system (1-2 users) for 
Informix k $3,800 for Unix and $995 for OS/2. Run-time prices are $1,540 for UNIX and $295 
for OS/2. The costs for Informix include the 4GL compiler. 


Object-Oriented Programming (OOP). Graphical user interface are composed of objects 
such as windows, buttons, etc. and are, therefore, a natural application for OOP techniques. 
Programming using objects allows the developer to work with objects such as windows 
while hiding the detaik of how the object works. This saves development and maintenance 
costs and facilitates portability between operating environments. The detaik of how an object 
works are buried in the definition of an object, not in the application that uses the object. 

Only the object library needs to be changed when porting between environments. The two 
OOP languages suited for this project are Smalltalk and C++. Smalltalk k a pure OOP 
environment while C++ k a hybrid language conskting of OOP extensions to C. Smalltalk V 
from Digitalk k available for Macintosh, Windows, and OS/2 environments. A UNIX 
version k planned. While Smalltalk k the best choice based on technical reasons, C++ k 
better for this project given the preference for C expressed by the peer review panel. The best 
C++ environment that k available on both UNIX and OS/2 k Glockenspiel C++ with the 
CommonView2 class library. The cost k $900 for OS/2 and $5,500 for IBM RS/ 6000 series 
machines. 
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Other Tools. There are several other tools available that generate executable files that require 
runtime licences. These tools are not discussed given the preference for C expressed by the 

peer review panel. 


6.3 RECOMMENDATIONS 

MTI recommends that we persue a two track development plan: 


■ Develop the scientific codes on a UNIX platform used by a majority of the 

large aerospace companies. IBM ADC can be used if eventual migration to 
the new IBM-Apple operating system is desired. 


Continue development of Industrial codes on an OS/2 platform with 
eventual migration to the new IBM-Apple operating system. 

Integrate the two systems using network communications. The integration 
task will be quite simple once DCE protocols are available for UNDC and 

OS/2. 


Object oriented tools such as Glockenspiel C++ and CommonView 2 class libraries wWchare 
portable between OS/2 and UNIX should be used to develop all the user interface code. Hus 
approach wifi leave open the option of using only one development platform and then 
porting to other platforms by recompiling the code on those platforms. In that case, 
recommends that OS/2 be the development platform given the reduced cost of developmen 
tools and lower development costs because of a simpler operating environment. 


ANSI standard FORTRAN77 should be used to develop the analytical codes. The codes 
should be designed such that they can also be used without the graphical user interfaces 
developed for the KBS. This will permit their distribution to customers using hardware and 
software platforms not supported by the KBS. 


Use a portable shell such as NEXPERT to develop the larger expert systems. Smaller 
embedded systems can be developed using CUPS or other C libraries. 


6.4 IMPLEMENTATION OF THE SEAL ANALYSIS KBS 

The work done during this period focused on the development of a graphical user mterface 
for some of the industrial seal codes and the executive program. OS/2 Presentation Manager 
(PM) and other system facilities were used to provide an interface with windows, drop-down 
menus, context sensitive help, dialog boxes for program input, and interactive graphics to 
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reduce the amount of numeric input where it was feasible to do so. The user interface was 
implemented for the following codes: 


■ The executive program that is used to control access to all the analysis 
codes and provide utility services such as printing and browsing text files. 

■ Spiral Groove Gas Seals Analysis (SPIRALG) 

■ Cylindrical Gas Seals Analysis (GCYL) 

■ Incompressible Cylindrical Seals Analysis (ICYL) 

■ Fluid Properties Calculations (FLUIDPROPS) 

■ Spiral Groove Face Seals Optimization Program (SPERALP) 


6.4.1 Description of User Interface 

The opening screen for the executive program is shown in Figure 6-2. Each program has its 
own button displaying the icon for the program. The user only has to click on the program 
button with a mouse to start a program. User options are selected from drop-down menus 
accessed from the action bar using either a mouse or a keyboard. For example, the FILE menu 
has options for printing and browsing text files such as output files created by the analysis 
programs. Figure 6-3 shows the output file from an analysis program being viewed in the 
browse window. Figure 6-4 shows the file selection screen for selecting input and output 
files. Plotting capability is currently provided using existing PC-DOS programs. A fully 
integrated OS/2 capability will be added in the future. 

Multitasking features of OS/2 are used to allow the user to have several codes running at the 
same time. The number of codes active at any time is limited only by memory available on 
the computer. Figure 6-3 shows die executive program and a Cylindrical Gas Seals analysis 
code (GCYL) active at the same time. The analysis code has been reduced to an icon (the 
same icon that is used in the button) to reduce screen clutter. Within each program, input 
and analysis are run as separate processes. This allows the user to start analyzing a data set 
and then prepare the input for the next analysis while the current analysis is in progress. 


Names of menu items are kept consistent between programs. For example, seal analysis 
variables are categorized according to function such as defining file scope of the analysis, 
specifying seal geometry, operating conditions and lubricant properties, etc. These functional 
groups are the same for most types of seals but the variables in each group change 
depending on the type of seal. The INPUT menu shown in Figure 6-5 lists the functional 
groups applicable for GCYL. Selecting a group from the list opens up a dialog box for 
entering values for variables in the group. All industrial codes have the same input menu list 
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FIGURE 6-2 Executive Program Main Window 
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FIGURE 6-3 Browsing Output File 


6-15 


V91-241 





FIGURE 6-4 File Selection Screen Input and Output 
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FIGURE 6-5 Standardized Input Menu Items 
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but the contents of the dialog boxes change depending on the seal type being analyzed. 
Examples of dialog boxes for the Analysis Options and Grid Definition menu items are 
shown in Figures 6-6 and 6-7, respectively. This consistency allows the user to quickly locate 
the variables to be input in any of the codes. Within each functional group, the names of 
variables have been made consistent across programs. Variable and user interface 
consistency should reduce learning time and make the codes easy to use by reducing the 
volume of information the user has to leam to master the interface. The user is left free to 
concentrate on the technical content of the codes. 


User input for analysis programs is done using dialog boxes containing entry fields for 
numeric data, radio buttons for selecting mutually exclusive options, and check boxes for 
selecting other optional features. The choices are presented in simple language avoiding 
computer jargon. Default values are provided for all variables. The user can move between 
fields using the mouse or the keyboard. Figure 6-8 shows a data entry screen from a code to 
calculate fluid properties. The input options are restricted to admissible values. For example, 
when the user selects fluid property calculations for specified temperatures and density, only 
the temperature and density entry field are displayed. The values entered in the fields are 
checked against acceptable limits. The limits are dynamic, and may change depending on the 
values of related variables. When the user switches the types of units used in the input, the 
values displayed in the entry fields and the unit labels are changed to reflect the choice of the 
user. The output is displayed using the same units as the input. The input values are saved 
by clicking on the ACCEPT button. Clicking on the DISCARD button discards the changes. 


Interactive graphics capability is provided where needed to reduce the amount of numeric 
input and to make the input more intuitive. For example, seal pads in padded seals analyzed 
by the GCYL code have several features on them such as recesses, Rayleigh steps, and fluid 
sources. In the original program, the user had to input the grid coordinates for the location 
and extent of all these features. An interactive capability is provided to enable the user to lay 
out the features on the grid using the mouse. Figure 6-9 shows a seal pad with a Rayleigh 
step and several constant pressure points. The user can add or delete the features shown in 
the features palette by using the mouse. The user first selects a feature from the palette by 
clicking on the appropriate radio button. The mouse pointer changes its shape to reflect user 
selection. The pointer is then moved to the grid window and the feature is placed on the grid 
by clicking the mouse at the appropriate grid location. Grid coordinates of the mouse pointer 
are displayed in the lower left comer of the grid display window. If additional information 
such as step height for a Rayleigh step is required for a given feature, entry fields are 
displayed above the features palette. The user interface code handles the details of generating 
the correct input statements for the GCYL code. 


Help is available at any time through the HELP menu or by pressing the FI function key. The 
FI key help is context sensitive. For example, if the user is entering data in an entry field and 
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Groove Gas Seals Program. 
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FIGURE 6-8 Data Input Methods 
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FIGURE 6-9 Interactive Specification of Seal Features 







presses the FI key, the help information for that field is displayed. Figure 6-10 shows the help 
window that pops up when the FI key is pressed while entering data in the "Groove Angle" 
field in GCYL. Figures from manuals are included in the help system. Once the help window 
is displayed, the user is free to browse through any portion of the help system for that code 
and has access to all the help utilities such as searching, printing, etc. provided by the 
Information Presentation Facility (IPF) in OS/2. Hypertext links are used as needed to 
provide explanations for technical terms used in the help information. The help for each 
program includes the following information: 

■ The purpose of the program 

■ Its capabilities and limitations 

■ References for additional information 

■ Code validation. 

■ Description of input and output parameters 

■ Examples describing the problem and showing typical input and output 
data sets 

■ Description of procedures for the user interface 


6.4.2 Current Status of KBS Components 

This section describes the status of the various components of the seal analysis KBS. 


6A.2.1 Executive Program 

The structure of the executive program is shown in Figure 6-11. The main program is 
designed to use separate threads for utility functions such as printing files and plotting data. 
This allows the user to have access to other program functions while these functions are 
being performed. When the user clicks on a program button, the analysis program is 
launched as a separate process. The button is disabled to prevent the use of multiple 
instances of the same program. The utility functions are available through the File menu. 
Utilities to print and browse output files have been implemented. Plotting capability will be 
added later. 


Help is available from the Help menu or by using the FI function key. The help information 
needs to be updated. 
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FIGURE 6-11: STRUCTURE OF THE EXECUTIVE PROGRAM 
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The name of the printer designated as the default printer during OS/2 setup is displayed in 
the main window. If the printer setup is changed, the display can be updated using the 
Printer Setup option in the Setup menu. Printer interface was designed for compatibility with 
OS/2 version 1.2. It needs to be updated to version 1.3 to provide additional capabilities that 
will make printing more flexible. The code to do so is already available but needs to be 
incorporated into the program. 


When the user selects the Print... option in the File menu, a file selection dialog box pops up 
to select the name of the file to be printed. Clicking on the Cancel button in the dialog box 
will cancel the printing procedure. After a file has been selected, a dialog box listing all 
available printer fonts is displayed to select the font to be used for printing. 


When the user selects the Browse... option in the File menu, a file selection dialog box pops 
up to select the name of the file to be browsed. Clicking on the Cancel button in the dialog 
box will cancel the browsing procedure. After a file has been selected, it is displayed in a 
separate window. The Fonts! menu item in file browse window lets the user select any of the 
available screen fonts. The window is closed using the system menu bar in the browse 
window. Only one browse window is allowed. 


The font support for browsing and printing needs to be improved. File opening and saving 
dialog boxes also need improvement. This work was postponed until OS/2 version 2 
becomes available because these dialog boxes have been standardized in that version. These 
changes are quite simple and will be implemented for all the codes. 


6A.2.2 Spiral Groove Gas Seals Analysis (SPIRALG) 

This code is complete. Enhancements made to the analysis codes after the user interface was 
developed may require some changes to the code. 

6.4.2.3 Cylindrical Gas Seals Analysis (GCYL) 

The help information in the code needs to be expanded. 


The user interface needs to support input for variable grids and to build in checks in the 
interactive seal layout portion of the code to prevent the user from specifying invalid seal 
configurations. The code was structured to support such checks but they need to be 
implemented. 
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6.4.2.4 Incompressible Cylindrical Gas Seals Analysis (ICYL) 


ICYL user interface was implemented using Toolbook. It is an alternate way of designing 
user interfaces for OS/2 applications. The menu options are displayed in the main window 
as shown in Figure 6-12 and are selected using a mouse or typing the number assigned to the 
item. Cascading menus shown in Figure 6-13 are used to display options available for each 
of the main menu items in a manner similar to drop down menus used in the presentation 
manager. An Input screen for menu item 5a is shown in Figure 6-14. Input elements such as 
radio buttons and entry fields are the same as for presentation manager applications. The 
user can go directly to any input screen using buttons in the upper right comer of the screen. 
Our work revealed some important shortcomings in Toolbook when used for engineering 
applications. The major difficulty was in the handling of large arrays. These can, however, be 
overcome by developing C functions which Toolbook can call. 


If NASA decides to stay with OS/2 for the CFD contract, file ICYL interface will be 
implemented using the presentation manager format because several prospective users at the 
last workshop at NASA wanted to stay with C. 

6.42.5 Fluid Properties Calculations (FLUIDPROPS) 

This code was obtained from NASA. The user interface is complete. 


Help information needs to be added to the program. 


The analytical portion of the program received from NASA is prone to crashes and was not 
changed in any way. Error trapping needs to be improved to facilitate graceful recovery from 
errors. 

6.42.6 Spiral Groove Face Seals Optimization Program (SPIRALP) 

The user interface is complete. 


Help information needs to be added to the program. 


The analytical portion of the program is prone to crashes and was not changed in any way. 
Error trapping needs to be improved to facilitate graceful recovery from errors. 
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FIGURE 6-12 User Interface for ICYL designed 
using Toolbook. 
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FIGURE 6-13 Cascading Menus in the ICYL Interface 
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FIGURE 6-14 Input Screen in the ICYL Interface. 
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6.5 FUTURE PLANS 


Additional analysis codes are currently being developed. The development of the user 
interface has been postponed pending final selection of an operating system. The decision is 
expected in October 1991. The development of expert system components will begin when 
tiie first scientific code is available in 1992. 
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